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bustion,  lean  blow-out,  and  thermo-acoustic  instabilities.  For  the  cases  of  stationary  com¬ 
bustion  and  lean  blow-out,  an  improved  version  of  the  Linear  Eddy  Model  approach  is  used, 
while  in  the  case  of  thermo-acoustic  instabilities,  the  effect  of  boundary  conditions  on  the 
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set  of  resolved  scale  Large  Eddy  Simulation  equations  for  continuity,  momentum,  energy,  and 
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is  that  the  Eulerian  species  concentration  equations  are  solved  at  the  resolved  scale  and  the 
Linear  Eddy  Model  is  strictly  used  to  close  the  species  production  term. 

In  this  work,  the  improved  Linear  Eddy  Model  approach  is  applied  to  predict  the  flame 
properties  inside  the  Volvo  rig  and  it  is  shown  to  over-predict  the  flame  temperature  and 
normalized  velocity  when  compared  to  experimental  data  using  a  premixed  single  step  global 
propane  reaction  with  an  equivalence  ratio  of  0.65.  The  model  is  also  applied  to  predict 
lean  blow-out  and  is  shown  to  predict  a  stable  flame  at  an  equivalence  ratio  of  0.5  when 
experiments  achieve  flame  extinction  at  an  equivalence  ratio  of  0.55.  The  improved  Linear 
Eddy  Model  is,  however,  shown  to  be  closer  to  experimental  data  than  a  comparable  reactive 
flow  simulation  that  uses  laminar  closure  of  the  species  source  terms. 
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CHAPTER  1 


Introduction 

In  air-breathing  propulsion  systems  it  is  vital  to  maintain  a  sustained  combustion  process 
in  fuel-air  flows  that  travel  at  speeds  exceeding  the  normal  burning  velocity  [1].  A  device 
that  has  been  shown  to  stabilize  and  maintain  a  flame  in  these  flows  is  a  bluff-body.  A  bluff- 
body  is  an  object  immersed  in  a  fluid  flow  that  causes  the  flow  to  travel  around  it,  creating 
highly  separated  flow  downstream  of  the  object  which  can  form  a  closed  recirculation  region. 
Bluff-bodies  have  been  commonly  used  in  premixed  combustion  applications  because  the 
recirculation  region  behind  the  object  serves  to  capture  heat  from  the  reactive  process  that 
continually  ignites  the  incoming  mixture  of  fuel  and  oxidizer  therefore  stabilizing  the  flame. 
It  is  important  to  understand  the  flow  fields  because  turbulent  bluff-body  stabilized  flames 
are  used  in  many  aerospace  propulsion  systems  including  gas  turbine  engines,  afterburners, 
ramjets,  and  scramjets. 

For  over  a  decade,  airbreathing  engine  development  has  emphasized  improving  efficiency 
and  reducing  emissions  as  the  driving  factors  behind  new  designs.  Recent  engines  have  been 
designed  to  run  leaner  in  an  effort  to  lessen  the  environmental  impact  of  emissions  such  as 
nitrogen  oxides  (NOx).  A  lean  engine  operates  with  an  O/F  ratio  less  than  the  stoichiometric 
ratio,  and  produces  lower  temperature  flames  that  reduce  thermal  nitric  oxide  (NO).  Lean 
combustion  also  emits  fewer  unburned  hydrocarbons  (UHCs)  because  the  excess  air  generally 
results  in  complete  consumption  of  the  fuel  [2],  Unfortunately,  combustion  in  this  lean  regime 
is  difficult  to  maintain  due  to  the  flames  being  more  susceptible  to  instabilities  and  lean 
blow  off,  especially  for  bluff  body  stabilized  flames.  Since  fuel  efficient  bluff  body  turbulent 
combustion  is  and  will  continue  to  be  important  to  many  aerospace  vehicles,  there  is  a  need 
for  accurate  prediction  of  turbulent  combusting  flows  in  this  regime.  This  need  has  led  to  the 
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research  and  development  of  improved  computational  fluid  dynamics  (CFD)  methods  and 
turbulent  combustion  models.  Improvements  of  these  methods  can  not  only  have  significant 
monetary  savings  but  also  serve  as  an  important  tool  for  the  design  of  new  efficient  engines 
that  operate  safely. 

This  thesis  will  review  several  turbulent  combustion  models  and  their  predictive  capa¬ 
bility  with  regards  to  bluff  body  stabilized  flames.  A  new  model  is  proposed  and  integrated 
into  an  established  CFD  code.  The  combined  solver  is  then  applied  to  a  bluff-body  stabi¬ 
lized  flame  problem  with  comparison  to  extensive  experimental  results  in  order  demonstrate 
its  viability  for  the  prediction  of  flow  properties  and  phenomena  such  as  lean  blow-out. 
Additionally,  a  study  is  also  conducted  on  how  boundary  conditions  and  experimental  geom¬ 
etry  affect  reacting  and  non-reacting  predictions  of  the  thermo-acoustic  interactions  inside 
a  two-dimensional  duct. 

1.1  Turbulent  Combustion  Modeling 

A  complete  simulation  of  all  relevant  spatial  and  temporal  scales  is  called  direct  numerical 
simulation  (DNS).  In  a  DNS  there  is  no  need  for  a  turbulence  or  turbulent  combustion  model 
[3].  All  spatial  scales  from  the  integral  scale  down  to  the  Kolmogorov  scale  are  captured 
on  the  computational  mesh.  Many  practical  aerospace  problems,  including  propulsion,  are 
associated  with  high  Reynolds  numbers  and  extremely  small  Kolmogorov  length  scales.  The 
number  of  grid  points  needed  is  excessive  and  therefore  computationally  very  expensive. 
DNS  also  requires  small  time  steps  that  adds  additional  expense.  As  a  result  of  this  required 
fine  resolution,  DNS  is  not  practical  for  most  relevant  aerospace  problems  and  is  therefore 
limited  in  its  application  to  academic  problems. 

Without  DNS,  some  type  of  model  is  required  to  capture  the  effects  of  the  turbulent  fine 
scales  that  are  not  being  directly  computed.  One  such  modeling  approach  uses  the  Reynolds 
Averaged  Navier-Stokes  (RANS)  equations,  which  have  long  been  used  in  CFD  simulations 
to  predict  turbulent  flows.  In  RANS  all  turbulence  scales  are  modeled,  which  makes  it  an 
attractive  model  because  relatively  coarse  grids  can  be  used.  The  coarse  grids  and  minimal 
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cost  of  the  RANS  simulation  make  it  computationally  tractable.  For  some  applications  it 
has  been  shown  to  produce  reasonable  results,  especially  for  steady,  non-reacting  flows.  It 
has  also  been  shown  to  be  unreliable  in  flows  involving  typical  propulsion  flowheld  conditions 
such  as  those  with  adverse  pressure  gradients,  flow  separation,  recirculation,  and  unsteady 
flame  dynamics  [4].  Underlying  reasons  for  the  inaccuracies  are  a  result  of  the  assumptions 
inherent  in  deriving  the  RANS  equations.  The  RANS  equations  are  derived  by  decomposing 
the  fluid  properties  (i.e.  pressure,  density,  velocity  and  species)  in  the  Navier-Stokes  equa¬ 
tions  into  mean  and  a  fluctuating  component.  The  resulting  equations  are  time  averaged  and 
contain  additional  terms  which  require  models  to  close  them  in  order  to  solve  the  equations. 
Many  different  models  such  as  the  k  —  e  and  the  k  —  to  [5]  models  have  been  proposed.  Each 
model  has  shown  success  for  specific  applications  but  no  model  is  universally  applicable,  gen¬ 
erally  because  they  employ  adjustable  constants  that  are  configuration  specific.  Additional 
modeling  is  required  for  reacting  flows  when  the  species  equations  are  solved  in  combination 
with  the  momentum  and  the  energy  equations.  It  is  clear  that  the  characteristics  of  the 
model  itself  becomes  vitally  important  to  the  accuracy  of  the  solution.  In  combustion,  the 
heat  is  released  at  the  small  scales  and  affects  the  larger  scales  so  it  is  important  to  cap¬ 
ture  the  interaction  between  these  scales.  In  RANS  models  there  is  no  distinction  between 
the  large  and  small  scales  which  means  that  the  modeled  terms  in  the  RANS  equations 
like  the  velocity-species  and  temperature-species  interactions  cannot  accurately  capture  the 
interaction  between  scales  [6]. 

A  natural  middle-ground  between  these  two  modeling  approaches  is  to  accurately  resolve 
a  portion  of  the  turbulence  and  model  only  the  smallest  scales.  The  large  eddy  simulation 
(LES)  framework  was  developed  with  this  in  mind.  LES  is  designed  so  that  it  resolves  the 
large  scale  turbulent  structures  and  models  the  small  scale  structures.  In  order  to  obtain  the 
scale  separation  in  a  reactive  flow,  the  governing  equations  and  species  equations  are  Favre- 
filtered[7]  which  then  allows  the  turbulent-kinetic-energy  containing  eddies  to  be  resolved 
and  the  smaller  scales  to  be  modeled.  LES  has  an  advantage  over  RANS  because  it  can 
represent  the  interaction  between  the  resolved  scale  and  the  unresolved  (modeled)  scale  [6]. 
The  filtering  operation,  for  LES,  like  in  Reynolds  averaging,  introduces  terms  that  are  not 
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closed  by  the  equations  and  therefore  need  turbulence  and  especially  turbulence  combustion 
models.  The  latter  refers  to  the  closure  of  the  species  source  term  in  order  to  accommodate 
the  influence  of  the  sub-grid  turbulent  fluctuations  on  the  combustion. 

There  are  many  turbulent  combustion  models  for  LES  that  have  been  examined  in  the 
literature.  Flamelet  [8,  9],  linear  eddy  modeling  (LEM)  [10,  11,  12,  13,  14,  15,  16],  and  the 
transported-probability  density  function  (PDF)  [17,  18,  19,  20]  are  three  of  the  most  widely 
used  models  in  the  literature  and  have  shown  success  in  modeling  different  phenomena. 
A  brief  description  of  each  of  these  sub-grid  models  will  be  discussed  to  identify  the  key 
underlying  assumptions  of  each  model. 

1.1.1  Steady  Laminar  Flamelet  Model 

The  steady  laminar  flamelet  model  employs  the  notion  that  a  turbulent  flame  front  can  be 
modeled  using  a  locally  one-dimensional  strained  flamelet  structure  inside  a  turbulent  flow. 
Flamelets  replace  the  species  and  energy  conservation  equations  with  a  single  conserved 
scalar  equation,  typically  in  terms  of  mixture  fraction.  Key  to  this  model  is  the  assumption 
that  the  chemical  kinetics  time-scales  are  smaller  than  the  smallest  turbulent  time  scales 
[8].  One  major  advantage  of  flamelets  is  that  the  conserved  scalar  equation  replaces  a  large 
number  of  species  transport  equations  and  is  therefore  less  computationally  expensive  to 
solve  compared  with  other  models.  Instead,  the  flamelet  approach  relies  on  a  pre-computed 
flamelet  library [9]  in  order  to  determine  local  species  mass  fractions,  thermodynamic  and 
transport  properties,  temperature,  and  source  terms. 

The  flamelet  model  involves  numerous  assumptions  which  limit  its  applicability.  Two 
key  assumptions  are  required  in  order  to  reduce  the  species  and  energy  equations  to  the 
conserved  scalar  form,  which  are,  a  constant  pressure  and  a  unity  Lewis  number,  where 
the  Lewis  number  is  the  ratio  of  thermal  diffusion  to  species  diffusion.  Flamelet  models 
are  therefore  only  valid  in  a  low  Mach  number  regime.  This  approximation  further  means 
that  the  thermal  diffusivity  and  all  species  mass  diffusivities  are  equal,  thereby  removing 
differential  mass  diffusion  which  is  important  for  many  reactive  flows.  Viscous  dissipation  is 
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also  explicitly  removed  from  the  conservation  equations  in  the  flamelet  modeling.  As  a  result 
of  these  assumptions,  unsteady  reactive  phenomena  and  coupling  to  acoustic  speeds  will  not 
be  correctly  predicted  by  the  model.  An  accurate  representation  of  acoustic  speeds  is  essential 
for  the  prediction  of  thermo-acoustic  instabilities  which  are  an  important  consideration  in 
the  design  and  evaluation  of  bluff-body  stabilized  propulsion  systems. 

The  flamelet  procedure  starts  with  a  pre-processing  step  where  the  flamelet  equations  are 
solved  and  tabulated  in  terms  of  mixture  fraction  and  local  scalar  dissipation  rate  [8].  Solving 
these  equations  and  tabulating  the  data  allow  for  rapid  access  to  the  value  for  species  mass 
fraction,  temperature,  and  thermodynamic  and  transport  properties.  The  flamelet  solutions 
are  based  on  a  canonical  configuration  like  a  counterflow  diffusion  flame  or  a  premixed 
flame[8].  Next,  the  source  term-free,  conserved  scalar  equation  is  solved  using  the  local  flow 
conditions  to  obtain  a  value  for  the  conserved  scalar,  which  is  a  parameter  constructed  from 
combinations  of  species  mass  fractions  and  temperature.  Using  the  calculated  conserved 
scalar  value  and  the  scalar  dissipation  rate  which  is  proportional  to  the  local  rate  of  strain, 
the  relevant  scalars  are  obtained  from  the  tabulated  data.  After  a  given  scalar  is  determined, 
it  then  needs  to  be  filtered  to  obtain  the  mass  fraction  variables.  In  order  to  perform  this  step, 
a  f3  joint  probability  density  function  (PDF)  function  for  the  conserved  scalar  is  typically 
used  [8].  The  reactive  scalar  and  conserved  scalar  are  then  integrated  to  obtained  a  filtered 
reactive  scalar  quantity. 

The  f3  PDF  is  a  further  limitation  of  the  flamelet  method  because  it  is  a  simplification 
compared  to  the  distribution  functions  of  other  turbulent  combustion  closure  methods  such  as 
LEM  and  transported  PDF.  Additionally,  there  is  some  confusion  about  flamelets  in  regions 
of  flow  that  have  both  premixed  and  non-premixed  combustion  because  the  correct  canonical 
configuration  to  use  in  each  domain  is  unclear.  Finally,  the  tabulated  data  are  computed 
based  upon  the  inlet  conditions  of  the  reactants  and  may  not  be  truly  representative  of  the 
local  flow  conditions. 
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1.1.2  Transported  Probability  Density  Function  Model 


The  transported  probability  density  function  method  models  the  non-linear  effects  in  turbu¬ 
lent  combusting  flows  using  a  one-point,  one-time  joint  PDF  for  a  given  set  of  flow  variables 
[17,  18,  19,  20].  The  PDF  method  was  originally  developed  to  be  used  as  a  RANS  closure 
and  more  recently  is  being  applied  to  LES  [21],  Currently,  the  most  advanced  PDF  methods 
incorporate  the  velocity,  the  turbulence  frequency,  and  the  composition  as  the  flow  variables. 
A  main  advantage  of  this  method  is  the  chemical  source  term  is  represented  exactly  in  the 
species  equations,  unlike  in  other  turbulent  combustion  models,  where  the  source  terms  need 
to  be  closed  by  introducing  other  model  assumptions.  These  source  terms  are  not  dependent 
on  spatial  or  temporal  correlations  and,  therefore,  they  are  ideal  to  be  represented  by  a 
single-point,  one-time  PDF  approach. 

The  PDF  equation  is  a  transport  equation  that  is  applicable  to  premixed  and  partially 
premixed  flames.  The  acceleration,  the  velocity  and  the  composition  (chemical  source  term) 
in  the  equation  are  solved  using  local  fluid  values.  The  viscous  effects,  the  fluctuating  pressure 
gradients,  and  the  molecular  diffusion  require  models  to  close  the  PDF  equation.  The  current 
area  of  research  in  PDF  based  methods  is  to  develop  improved  models  to  capture  the  mixing 
and  diffusion  [22,  23].  The  transported  PDF  method  is  modeled  using  a  low  Mach  number 
or  constant  pressure  approximation  that  limits  its  applicability  to  many  aerospace  relevant 
flows. 

The  PDF  equation  is  solved  using  Lagrangian  methods  due  to  the  large  number  of  inde¬ 
pendent  variables  involved  in  these  equations.  Particle-based  method  such  as  Monte  Carlo 
methods  are  used  because  finite  element  and  volume  methods  are  much  more  computation¬ 
ally  expensive  [17].  The  most  significant  disadvantage  of  the  PDF  method  is  that  the  large 
number  of  independent  variables  involved  in  solving  the  equations  make  this  method  com¬ 
putationally  expensive.  In  addition,  obtaining  a  proper  mixing  model  for  reacting  flows  has 
been  a  limitation  of  this  sub-grid  model. 
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1.1.3  Linear  Eddy  Model  (LEM) 


The  linear  eddy  modeling  approach  was  developed  by  Alan  Kerstein  et  al.  [10,  11]  to  predict 
scalar  mixing  in  turbulent  flows  by  stochastically  representing  the  effect  of  turbulence  on 
a  scalar  profile  at  all  relevant  physical  length  scales.  LEM  was  first  tested  as  a  turbulence 
model  and  produced  accurate  results;  it  was  later  adapted  to  model  turbulent  combustion 
[12,  13].  The  LEM  approach  is  a  multi-scale  model  that  embeds  a  one-dimensional  grid 
in  the  sub-grid  of  the  multi-dimensional  CFD  mesh.  The  LEM  grid  is  used  to  solve  the 
conservation  equations  for  the  species  and  the  low  Mach  number  form  of  the  energy  equation. 
The  orientation  of  the  LEM  line  is  such  that  it  is  parallel  to  the  maximum  species  gradient 
within  each  LES  cell.  The  three  physical  phenomena  described  by  the  LEM  equations  are 
turbulent  stirring/advection,  molecular  diffusion,  and  combustion.  An  attractive  feature  of 
the  LEM  approach  is  that  the  full  species  equations  are  solved  and  not  approximated  as  in 
other  sub-grid  models,  although  the  solution  is  obtained  in  one-dimension. 

A  unique  feature  of  LEM  is  that  the  turbulent  stirring  is  modeled  stochastically  using  a 
rearrangement  process  called  a  triplet  map  [10].  The  triplet  map  simulates  the  effect  that  the 
inertial  range  turbulence  has  on  a  scalar  profile  [15].  The  map  shrinks  the  one-dimensional 
scalar  profile  into  a  third  of  the  original  length  and  then  fills  the  interval  with  two  additional 
compressed  copies  of  the  profile.  The  mapping  process  mimics  the  effect  of  a  vortex  crossing 
while  maintaining  complete  species  conservation.  The  frequency  and  the  location  of  the 
stirring  event  are  chosen  stochastically  using  probability  distribution  functions  [11].  The 
molecular  diffusion  in  the  LEM  is  accounted  for  using  using  Fick's  law  with  the  Hirschfelder 
and  Curtiss  approximation,  wherein  a  correction  velocity  is  added  to  the  molecular  diffusion 
to  ensure  conservation  of  mass  for  all  species  [7]. 

LEM  is  implemented  as  the  sub-grid  model  for  LES  to  form  the  LEM-LES  method. 
LEM-LES  solves  the  filtered  form  of  the  continuity,  momentum  and  energy  equations  at  the 
resolved  scale  and  solves  LEM  at  the  unresolved  scale.  It  solves  separate  energy  equations  at 
the  resolved  and  unresolved  scales.  Solving  two  energy  equations  can  produce  an  inconsistent 
temperature  between  the  two  scales.  In  the  sub-grid,  a  low  Mach  form  of  the  energy  equation 
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implies  that  the  pressure  is  constant  in  each  LEM  cell.  The  temperature  calculated  at  the 
sub-grid  level  is  erased  after  each  physical  time  step,  and  serves  only  to  integrate  the  species 
transport  equation  along  the  LEM  line. 

LEM-LES  does  not  solve  the  Eulerian  form  of  the  species  transport  equations  in  the 
resolved  scale.  Rather,  it  uses  a  process  called  Lagrangian  splicing  to  account  for  the  species 
transport  between  LES  cells  [14].  Splicing  divides  the  mass  flux  of  each  LEM  line  into 
individual  inflow  and  outflow  contributions  so  the  quantity  of  mass  entering  and  leaving 
each  LEM  line  is  tracked.  The  mass  transfer  is  accomplished  by  using  the  calculated  outflux 
to  remove  the  corresponding  number  of  LEM  cells  from  the  output  end  of  the  LEM  line. 
Next,  the  mass  influx  is  used  to  add  the  outflux  LEM  cells  from  other  LEM  lines  onto 
the  input  side  of  the  LEM  line.  After  the  mass  has  been  transferred  between  LES  cells, 
the  number  of  LEM  cells  within  the  LEM  domain  has  changed  as  a  results  of  the  splicing. 
Finally,  the  LEM  domain  is  re-gridded  in  order  to  maintain  a  uniform  LEM  mesh  [24], 

It  is  important  to  stress  that  in  LEM-LES  the  Eulerian  species  equations  are  not  solved  at 
the  resolved  scale  but  instead  the  species  transport  is  accounted  for  solely  in  LEM  model.  The 
Lagrangian  splicing  process  has  an  adherent  randomness  due  to  the  arbitrary  distribution 
of  mass  between  LES  cells.  This  random  distribution  of  mass  occurs  because  the  one¬ 
dimensional  LEM  line  is  interacting  in  a  multidimensional  space.  Additionally,  the  large  scale 
advection  process  only  accounts  for  advection  of  mass  and  does  not  incorporate  the  large 
scale  diffusion  that  takes  place  between  LES  cells  which  means  that  the  species  equations 
do  not  converge  to  the  correct  DNS  solution  in  the  limit  of  very  fine  grids. 

1.2  Bluff-Body  Stabilization  of  Flames 

As  noted  earlier,  flames  stabilized  by  bluff-bodies  are  used  in  many  aerospace  propulsion 
systems.  As  a  result  there  have  been  numerous  computational  studies  to  assess  the  predictive 
capability  of  CFD  tools  to  predict  key  flow  physics.  The  difficulties  in  these  problems  arise 
due  to  the  complex  vortex  dynamics  that  are  seen  during  stable  and  unstable  combustion 
[25,  26].  Additionally  the  challenge  of  predicting  phenomena  such  as  lean  blow-out  and 


thermoacoustics  provides  added  difficulty  to  state-of-the-art  CFD  codes  because  they  are 
inherently  unsteady  phenomena  that  have  a  strong  dependence  on  local  flow  conditions 
and  can  often  not  be  captured  by  steady  state  simulations.  The  challenges  involved  in 
computationally  predicting  the  physical  flow  parameters  arise  because  the  CFD  solution  is 
greatly  influenced  by  the  numerical  scheme  and  the  turbulent  combustion  closure.  Efforts  to 
improve  CFD  predictions  of  bluff-body  stabilized  flames  have  been  tested  using  a  number  of 
canonical  configurations  such  as  the  Volvo  rig  [27,  28,  29]  which  is  a  rectangular  duct  that 
contains  an  equilateral  triangle-shaped  bluff-body  flame  holder. 

The  effect  of  the  numerical  discretization  scheme  on  the  predictions  of  the  Volvo  rig 
flowheld  was  studied  by  Cocks  et  al.  [30],  and  the  results  showed  that  grid  dependency 
and  numerical  errors  can  greatly  impair  the  predictive  capabilities  of  the  LES.  The  results 
showed  that  despite  using  the  same  grid,  turbulence  model,  and  flow  parameters,  different 
numerical  schemes  (all  of  which  are  formally  second  order  accurate)  can  produce  remarkably 
different  time-dependent  results. 

The  three  main  phenomena  of  interest  in  bluff-body  stabilized  flames  are  stationary 
combustion,  lean  blow-out,  and  thermo- acoustic  instabilities.  The  accurate  prediction  of 
flow  properties  in  these  regimes  is  heavily  dependent  on  the  turbulent  combustion  model 
used  in  the  solver.  As  a  result,  there  has  been  much  focus  in  past  research  on  applying 
a  wide  range  of  turbulent  combustion  models  to  bluff-body  stabilized  flames  to  asses  their 
validity  in  predicting  flow  properties  in  these  three  regimes. 

Turbulent  combustion  modeling  studies  have  mainly  focused  on  prohle  comparisons  of  the 
root-mean-sqnare  values  of  the  fluctuating  temperature  and  the  velocity  in  stable  premixed 
flows.  One  of  the  earliest  studies  to  predict  flow  inside  the  Volvo  configuration  using  LES  was 
performed  by  Fureby  and  Moller  [31].  The  goal  of  this  study  was  to  investigate  the  effect  that 
sub-grid  models  have  on  the  predictive  capabilities  of  LES.  Their  results  indicated  that  LES 
is  capable  of  predicting  most  significant  flow  features,  including  unsteady  flow  structures. 

Flamelets  are  an  attractive  choice  for  combustion  modeling  due  to  the  low  computa¬ 
tional  cost  compared  to  many  other  models.  Fureby  [32]  applied  both  the  propagation-based 
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flamelet  model  and  the  filtered  flam  el  et,  model  to  premixed  bluff  body  stabilized  flames.  The 
results  showed  that  both  of  the  flamelet  models  produced  reasonable  comparison;  however, 
the  filtered  flamelet  model  was  more  sensitive  to  the  grid  resolution.  Additional  research  in 
flamelets  was  conducted  by  Baudoin  et  al.  [33]  in  which  the  G-equation  model  and  four  finite 
rate  chemistry  models  were  employed.  These  chemistry  models  include  the  partial  stirred 
reactor  model  (PaSr),  thickened  flame  model  (TFM),  eddy  dissipation  concept  (EDC),  and 
presumed  probability  density  function  (PPDF).  These  were  compared  using  the  Volvo  rig 
geometry.  The  results  showed  that  all  four  finite  rate  models  had  very  good  agreement  with 
each  other  and  experimental  data  in  the  temperature,  the  velocity,  and  the  species  concen¬ 
tration  despite  being  significantly  different  modeling  approaches.  On  the  other  hand,  the 
flamelet  G-equation  model  had  a  more  noticeable  deviation  with  experimental  data  compared 
to  the  finite  rate  models. 

Porumbel  and  Menon  compared  the  eddy  break-up  (EBU)  and  LEM  as  closure  conditions 
for  the  LES  resolved  scale  equations  [34],  The  simulations  were  run  with  both  reacting  and 
non-reacting  flow  using  the  Volvo  rig  geometry.  The  reacting  simulation  used  a  premixed 
propane-air  mixture  at  an  equivalence  ratio,  0,  of  0.65.  Both  models  predicted  Von  Karman 
vortex  shedding  and  other  flow  parameters  accurately  for  the  non-reacting  flow  case.  For 
the  reacting  flow  case,  both  models  predicted  the  time-averaged  velocity  with  good  accu¬ 
racy,  while  the  EBLI  sub-grid  model  failed  to  predict  the  turbulent  fluctuations  in  the  flow 
field.  Additionally,  the  EBU  under-predicts  the  centerline  temperature  due  to  insufficient 
prediction  of  turbulent  mixing  rate. 

Further  work  was  done  by  Jones  et  al.  [35]  using  a  sub-grid  joint  PDF  in  conjunction  with 
the  Eulerian  stochastic  field  method  as  the  turbulent  combustion  closure.  The  results  showed 
good  agreement  with  the  experimental  results  and  showed  that  LES  accurately  predicts 
premixed  combustion  fluid  properties  when  coupled  with  an  appropriate  sub-grid  model. 
From  analyzing  past  work  it  is  possible  to  conclude  that  accurately  predicting  the  mean  flow 
parameters  (temperature,  velocity,  species  concentration)  in  a  stable  premixed  turbulent 
combustion  problem  is  possible  with  a  wide  number  of  closures. 

Research  into  lean  blowout  prediction  in  premixed  bluff  body  stabilized  flames  was  per- 
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formed  by  Gokulakrishnan  et  al.  [36]  on  the  experimental  setup  by  Kiel  et  al.  [37]  at  the  Air 
Force  Research  Laboratory  (AFRL).  This  experiment  differs  from  the  Volvo  experiment  in 
that  the  bluff  body  is  a  V-gutter  instead  of  an  equilateral  triangle.  CFD  calculations  were 
performed  using  the  laminar  chemistry  model  (LCM)  [38]  and  EDO  [39]  models  to  determine 
the  equivalence  ratio  that  computationally  predicts  lean  blowout.  The  results  showed  that 
LES  coupled  with  EDC  showed  blow-out  at  an  equivalence  ratio  of  0.6,  which  matches  the 
observed  experimental  results.  The  laminar  chemistry  model  did  not  blowout  until  a  much 
lower  equivalence  ratio  of  0.45.  This  research  demonstrates  that  the  turbulence-chemistry 
interactions  play  a  key  role  in  predicting  the  wake  region  behind  the  flame  holder,  and 
therefore  flame  stability. 

The  effect  of  boundary  conditions  and  domain  configuration  was  researched  by  Ma  et  al. 
[40]  on  the  Volvo  rig  to  analyze  the  prediction  of  thermo-acoustic  phenomena.  Two  simula¬ 
tions  were  run,  one  with  a  full  Volvo  rig  domain  and  a  constant  pressure  outlet  boundary 
condition,  the  other  simulation  used  a  partial  domain  with  an  impedance  boundary  condi¬ 
tion  for  the  outlet.  The  simulations  were  performed  using  the  FLUENT  commercial  software 
package  with  the  flamelet  model  being  used  for  turbulent  combustion  closure.  The  focus  of 
the  work  was  to  compare  the  pressure  wave  amplitude  and  frequency  content  for  the  two  dif¬ 
ferent  outlet  boundary  conditions.  The  results  showed  that  both  cases  had  good  agreement 
in  terms  of  the  predicted  frequency  but  the  magnitude  of  the  oscillations  was  different  for 
some  equivalence  ratios.  The  results  also  showed  an  increase  in  pressure  oscillations  with  an 
increase  in  equivalence  ratio. 

1.3  Computational  Fluid  Dynamics  Solver 

The  CFD  code  used  in  the  present  research  is  the  General  Equation  and  Mesh  Solver,  or 
simply  GEMS.  GEMS  is  a  CFD  code  that  was  designed  by  Charles  Merkle’s  group  at  Purdue 
LIniversity  and  is  actively  being  used  by  Purdue  and  the  Air  Force  Research  Laboratory  to 
perform  both  reacting  and  non-reacting  CFD  calculations  in  rocket  and  turbine  engines 
[41,  42],  The  code  uses  a  finite  volume  method  to  solve  the  continuity,  the  Navier-Stokes, 
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and  the  energy  equations  with  second  order  accuracy  in  both  time  and  space.  In  addition, 
GEMS  is  able  to  handle  turbulence  and  combustion  by  solving  the  species  conservation 
equation  and  turbulence  model  equations.  For  the  turbulence  model  equations  GEMS  uses 
the  k  —  u  model  which  describes  the  turbulent  kinetic  energy  and  turbulent  dissipation  in 
the  flow,  respectively.  The  combustion  is  modeling  using  a  laminar  flame  rate  model  that  is 
a  function  of  the  mean  temperature,  the  pressure,  and  the  species  concentration  [42], 

GEMS  uses  a  combination  of  RANS  and  LES  called  Detached  Eddy  Simulation  (DES) 
[42]  to  capture  turbulent  flow  fields.  DES  uses  RANS  in  the  near- wall  boundary  layer  region, 
where  the  RANS  model  performs  efficiently  in  cells  with  high  aspect  ratios.  On  the  other 
hand,  LES  is  implemented  in  areas  away  from  the  boundary  layer  where  the  grid  resolution 
supports  the  turbulent  scales  [42],  thereby  capturing  the  unsteady  flame  dynamics  in  regions 
where  combustion  is  dominant. 

1.4  Thesis  Objectives 

The  goal  of  this  thesis  is  to  analyze  the  three  main  aspects  of  bluff-body  stabilized  flames; 
stationary  combustion,  lean  blow-out,  and  thermo-acoustic  instabilities.  In  the  cases  of 
stationary  combustion  and  lean  blow-out,  an  improved  version  of  the  LEM  model  is  used, 
while  in  the  case  of  thermo-acoustic  instabilities,  the  effects  of  boundary  conditions  on  the 
predictions  are  studied.  LEM  was  chosen  because  it  has  validity  for  all  flame  regimes  and 
is  a  multi-scale  approach  that  solves  the  full  set  of  species  equations  at  the  sub-grid  level, 
albeit  in  one- dimension. 

This  thesis  considers  an  improved  formulation  of  the  LEM  by  coupling  it  with  a  full-set  of 
resolved  scale  LES  equations  for  continuity,  momentum,  energy,  and  species  transport.  The 
main  difference  from  past  implementations  is  that  the  species  concentration  equations  are 
solved  at  the  resolved  scale,  which  includes  diffusion  and  advection.  In  the  past,  the  species 
concentration  calculations  have  been  completely  confined  to  the  sub-grid  LEM  model  and  the 
transport  of  species  between  LES  cells  only  accounted  for  advection.  This  implementation 
uses  LEM  to  obtain  a  closure  value  for  the  average  species  production  term  in  the  filtered 
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species  equations  at  the  resolved  scale. 

In  addition  to  improving  LEM,  this  thesis  analyzes  the  effect  of  boundary  conditions  and 
facility  geometry  for  studying  thermo-acoustics  in  bluff  body  stabilized  flames.  Many  CFD 
studies  have  been  performed  previously  using  the  Volvo  rig  geometry  but  they  have  typically 
not  modeled  the  farfield  downstream  of  the  duct  exit.  Instead,  a  boundary  condition  has 
been  set  at  the  duct  exit.  The  present  study  seeks  to  determine  if  physical  and  computational 
boundary  conditions  affect  the  solution  inside  the  duct  by  modeling  the  exhaust  region  for 
both  reacting  and  non-reacting  fluid  flows. 

The  chapters  for  the  thesis  are  organized  as  follows: 

•  Chapter  2  provides  a  detailed  description  of  the  LEM  approach.  First,  the  general 
equations  are  derived.  Next,  the  numerical  implementation  is  discussed.  Finally,  the 
assumptions  used  in  LEM  are  explicitly  stated  so  the  applicability  of  LEM  can  be 
established. 

•  Chapter  3  tests  the  predictions  of  the  LEM  model  for  a  single  sub-grid  cell.  This 
section  focuses  on  verifying  that  the  model  has  been  properly  implemented  and  that 
diffusion,  stirring,  and  combustion  processes  function  correctly. 

•  Chapter  4  uses  the  integrated  LEM-LES  solver  to  simulate  flow  in  the  Volvo  rig.  The 
results  for  temperature,  velocity,  and  lean  blow-out  equivalence  ratio  are  compared 
with  experimental  data. 

•  Chapter  5  analyzes  the  effect  of  boundary  conditions  and  facility  geometry  on  bluff 
body  stabilized  flames  using  the  geometry  of  an  experimental  rig  at  AFRL. 

•  Chapter  6  concludes  the  research  by  assessing  the  advantages  and  disadvantages  of  the 
improved  LEM.  Lastly,  a  proposal  for  future  LEM  improvements  is  made. 
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CHAPTER  2 


LEM  Approach 


The  filtered  species  production  term  in  the  LES  species  transport  equation  is  difficult  to 
calculate  in  turbulent  combustion  because  it  is  no  longer  a  function  of  P,  T,  and  Y  as  in 
laminar  combustion  and  therefore  requires  a  model  for  closure.  The  linear  eddy  model  is  a 
multi-scale  model  used  to  obtain  the  filtered  species  production  term.  The  implementation 
of  the  LEM  presented  in  this  work  differs  in  some  key  ways  from  other  implementations  in 
the  literature  and  is  therefore  described  in  detail  in  this  chapter. 


2.1  LES  Equations 


Large  eddy  simulations  (LES)  solve  the  filtered  conservation  equations.  The  Favre-filtered 
form  of  the  conservation  laws  for  mass,  momentum,  energy,  and  species  transport  equations 
are  [7]: 


dp  d  . 

at  +  ipUi)  = 0 

d  d  _ _  dp  d  _ _  __  __ 

Ft  ^  +  FFj  ^ pu jUi)  =  ~F^l  +  FF^fji~p  (UjUi  ~  UjUi ^ 


d_ 

Ft 

d_ 

Ft 


+ 


+ 


d 

dXj 

d 

dXn 


dp  d 

=  Ft +  FF-{ UiTij  ~  pi  +  ^ UiTij  ~ 

_  d 

—  +  "7 - 

OX, 


—  p  ( Ujho  —  Ujho 


pVk,jYi  H-  p  [ujYk  UjYk 


(2.1) 

(2.2) 

(2.3) 

(2.4) 


where  oji  is  the  filtered  species  production.  The  terms  with  a  tilde  are  mass-averaged  variables  and 
the  terms  with  the  bar  are  time  averaged.  The  Wilcox  k  —  u  turbulence  model  is  [5]: 
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These  are  the  equations  used  in  GEMS  (ref:  Section  1.3).  The  filtered  species  production  term  is 
closed  using  the  LEM  sub-grid  model  and  is  described  in  the  following  sections. 

2.2  Introduction  of  LEM 

LEM  implements  a  one-dimensional  form  of  the  species  mass  fraction  and  energy  conservation 
equations  to  describe  advection,  molecular  diffusion,  and  the  production  of  species  in  the  sub-grid 
of  an  LES  domain.  LEM  is  implemented  as  a  multi-scale  approach  by  solving  the  species  and 
energy  equations  at  both  the  resolved  and  unresolved  scales.  A  unique  feature  of  the  LEM  model  is 
that  the  sub-grid  equations  are  represented  by  a  one-dimensional  set  of  equations.  While  diffusion 
and  reaction  are  solved  by  resolving  down  to  the  Kolmogorov  scales,  sub-grid  advection  is  modeled 
using  a  mapping  process  that  simulates  turbulence  along  a  one-dimensional  domain  [10,  11]. 

The  effect  of  eddies  smaller  than  the  LES  filter  width,  A,  is  not  captured  in  the  resolved  scale. 
In  the  LEM  model  they  are  accounted  for  in  the  sub-grid  stirring  term,  Fsfir.  The  model  is  a 
rearrangement  process  called  a  triplet  map  that  simulates  the  effect  of  a  turbulent  eddy  passing 
through  a  one-dimensional  scalar  profile,  as  seen  in  Figure  2.1.  This  process  can  be  viewed  as 
an  instantaneous  event  that  compresses  the  scalar’s  spatial  profile  (Fig.  2.2a)  to  one-third  of  its 
original  size  (Fig.  2.2b).  The  compressed  profile  is  then  replicated  twice  to  fill  the  remaining 
domain  (Fig.  2. 2d).  To  conserve  mass  and  energy,  the  center  segment  is  inverted  (Fig.  2.2c).  The 
triplet  map  is  implemented  using  three  parameters:  the  local  eddy  size,  £,  the  frequency,  /,  and 
the  stirring  location,  so-  The  eddy  size  and  location  are  the  same  for  all  species  and  temperature 
profiles  during  a  single  time  step  so  that  the  model  remains  physically  consistent. 

2.3  Derivation  of  the  LEM  Equations 

In  LEM,  the  energy  and  species  equations  are  solved  in  the  sub-grid  on  a  one-dimensional  line  with 
a  spatial  resolution  commensurate  with  the  Kolmogorov  scale.  This  model  captures  turbulent- 
chemistry  interaction  by  modeling  the  molecular  diffusion,  advection,  and  the  production  of  species 
in  the  sub- grid.  The  traditional  equations  for  species  and  energy  conservation  are  modified  to 
enable  this  process.  The  details  of  these  equations  are  given  in  the  subsequent  subsections. 
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Figure  2.1:  Turbulent  eddy  acting  on  scalar  profile.  Adapted  from  [10]. 

2.3.1  Sub-Grid  Species  Equation 

The  quantities  of  interest  are  in  the  sub-grid,  therefore  it  is  appropriate  to  work  with  the  unfiltered 
conservation  equations.  The  unfiltered  species  equation  is 

+  {fiUiYk)  =  -£r,  {pV‘-kYk) +  Clk  (2'7) 

where  Vj  k  is  the  diffusion  velocity  and  ojk  is  the  species  production  term.  The  advection  velocity 
term,  Uj,  is  split  into  two  terms 

Uj  =  u  +  u"  (2-8) 

where  u  is  the  filtered  velocity  used  at  the  resolved  scale  and  neglected  from  the  sub-grid  equation. 
The  u"  term  represents  eddies  at  the  unresolved  scale  and  it  is  therefore  used  in  the  LEM  species 
sub-grid  equation 
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Equation  2.9  represents  the  sub-grid  species  equation  and  is  solved  along  the  LEM  line.  The 
diffusion  velocity  in  this  equation  is  modeled  using  Fick’s  Law  of  diffusion  with  the  Hirschfelder 
and  Curtiss  approximation  for  mass  diffusivity 

vi*y„  = 
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(2.10) 


(a)  Initial  profile. 


(c)  Reverse  of  compressed  profile. 


Figure  2.2:  The  triplet  map  process  is  shown 
pressed  to  one-third  of  the  original  length  (b), 
third  compressed  profile  is  added  to  complete 


(b)  Profile  compressed  to  1/3  of  length. 


(d)  Three  compressed  profiles. 


here.  The  original  scalar  profile  (a)  is  com- 
the  compressed  profile  is  mirrored  (c),  and  a 
die  mapping  process  (d). 
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This  equation  is  only  correct  for  binary  diffusion.  For  more  than  two  species  it  is  only  an  ap¬ 
proximation  which  can  lead  to  errors  in  mass  conservation.  To  eliminate  these  errors  a  correction 
velocity,  V?,  is  introduced.  The  final  form  of  the  species  equation  being  solved  in  LEM  is 
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where  FstiT  is  a  model  term  to  account  for  the  effect  of  the  sub-grid  eddies.  The  correction  velocity 
is  defined  as 


yc  _  Pk,M  d 

3  Wmix  dXj 


( WmixYk ) 


(2.13) 


for  Ns  species  comprising  the  mixture.  The  prime  notation  has  been  dropped  for  simplicity. 


2.3.2  Sub-Grid  Temperature  Equation 


The  unfiltered  energy  equation  with  no  body  forces  and  no  heat  addition  for  a  calorically  perfect 
gas  is 
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where  the  term  q%  is  the  enthalpy  flux.  As  before  the  velocity  can  be  decomposed  into  an  filtered 
value  (resolved)  and  a  fluctuating  component  (sub-grid).  The  resolved  scale  velocity  is  neglected 
from  the  sub-grid  equations  and  the  fluctuating  component  is  used. 

The  equation  is  written  in  terms  of  temperature  [7],  which  results  in  the  LEM  energy  equation 
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where  A  is  the  thermal  conductivity  and  wy  is  the  energy  addition  due  to  combustion  given  by, 
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where  hStk  is  the  sensible  enthalpy  and  hf  j.  is  the  heat  of  formation.  As  is  done  for  the  species 
equations,  the  advection  term  will  be  modeled  using  a  triplet  map  and  Fick’s  Law  for  diffusion  is 
implemented  for  the  diffusion  velocity.  The  final  form  of  the  energy  equation  is 
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where  TrjStir  is  a  model  term  to  account  for  the  effect  of  the  sub-grid  eddies. 
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2.4  Implementation  of  LEM 


LEM  embeds  a  one-dimensional  line  within  a  multi-dimensional  LES  grid.  The  LEM  line  should 
be  oriented  so  that  it  captures  the  important  flow  features  within  each  LES  cell.  In  other  LEM 
implementations  the  line  is  typically  oriented  to  align  with  the  maximum  species  mass  fraction 
gradient  [10].  For  this  work  the  LEM  line  is  instead  oriented  across  the  flame  front  which  can  be 
detected  by  examining  the  temperature  gradient,  as  described  in  section  2.4.2.  This  means  that 
the  orientation  of  the  line  will  be  different  in  each  LES  cell.  The  final  equations  are  then  written 
in  terms  of  the  parameterized  spatial  coordinate,  s,  and  are  given  as 
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2.4.1  LEM  Grid 


The  LEM  grid  is  a  one-dimensional  line  that  has  a  length  equal  to  the  LES  filter  width,  A.  It  is  im¬ 
portant  that  the  grid  have  sufficient  resolution  to  capture  turbulent  effects  down  to  the  Kolmogorov 
scale,  r/.  The  Kolmogorov  length  scale  is  defined  as  [43] 

7]  =  NriAR,e~^i/4  (2.20) 


where  Nv  is  an  empirical  constant  that  is  set  to  unity  in  this  work.  The  Reynolds  number  is  defined 
as 

ReA  =  ^  (2.21) 

where  v!  is  the  LES  sub-grid  velocity  field  and  is  defined  as 
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where  kS9S  is  the  sub-grid  kinetic  energy  obtained  from  the  resolved  scale.  To  resolve  the  Kol¬ 
mogorov  length  scales  the  grid  spacing,  As,  needs  to  be  smaller  than  ?/.  In  practice  the  grid  size 
must  be  three  times  smaller  than  the  Kolmogorov  scale  in  order  to  exactly  conserve  mass  during 
the  triplet  mapping  process.  Therefore,  the  number  of  LEM  cells  is  given  by 

Nlem  =  3  •  round  ^ ^  (2.23) 
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The  grid  spacing  is  determined  as: 


As 

A  schematic  of  the  LEM  grid  is  shown  below, 


A 

Nlem 


(2.24) 


LEM  node  LEM  ce|. 


- AS 

-L  =  A 

Figure  2.3:  LEM  grid  schematic. 


2.4.2  Initialization  of  LEM  solver 

In  this  research  the  LEM  line  is  aligned  with  the  resolved-scale  temperature  gradient  so  that  the 
flame  front  is  captured.  This  is  accomplished  by  using  the  temperature  gradient  to  determine 
the  normal  direction  of  the  LEM  line.  The  normal  vector  is  dotted  with  the  individual  species 
gradients  to  obtain  the  species  mass  fractions  along  the  line,  and  the  temperature  gradient  is  used 
to  determine  the  temperature  profile  along  the  line.  The  species  mass  fractions,  temperature, 
species  gradients,  and  temperature  gradient  are  obtained  from  the  resolved  scale  (i.e.  the  LES 
solution)  at  the  start  of  each  LES  time  step. 


The  magnitude  of  the  temperature  gradient  is  calculated  as  follows: 

|VT|  = 


IdT 2  dT' 


dT ‘ 


dx  +  dy  dz 


The  species  mass  fraction  gradients  along  the  LEM  line  are  determined  by: 


(2.25) 


VYk,llVT-  —  -VYk  (2.26) 

The  gradients  are  used  to  linearly  extrapolate  the  species  mass  fraction  and  temperature  profiles 
for  the  LEM  line  using  the  LES  value  at  the  center  of  the  line.  However,  the  extrapolation  process 
can  potentially  cause  the  species  mass  fraction  values  to  be  greater  than  unity  or  less  than  zero  at 
the  endpoints  of  the  LEM  cell.  This  situation  occurs  when  a  species  has  a  zero  mass  fraction  in 
one  cell  and  the  neighboring  cell  has  a  non-zero  value.  Any  gradient  in  this  situation  will  cause  half 
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of  the  profile  to  be  less  than  zero  because  the  line  is  centered  around  zero.  This  issue  is  mitigated 
by  automatically  setting  a  zero  gradient  for  any  species  that  has  a  zero  species  mass  fraction.  In 
doing  so,  the  total  species  mass  fraction  at  each  LEM  node  will  no  longer  add  up  to  unity,  except 
at  the  center  node.  It  is  important  in  the  initialization  process  that  the  total  species  over  the  LEM 
line  is  conserved,  meaning  that  each  nodal  sum  of  mass  fractions  is  unity.  In  order  to  achieve  mass 
unity  at  each  node  it  is  necessary  that  the  species  gradients  sum  to  zero  as  shown  below: 

Ns 

E^-  =  1  (2-27) 

k= 1 


NS  ATT 

k= 1 


ds 


=  0 


(2.28) 


In  order  to  meet  the  requirement  of  equation  2.28  the  gradients  need  to  be  adjusted  to  account  for 
any  species  gradient  that  was  set  to  zero  due  to  zero  mass  fraction.  The  largest  species  value  of 
the  remaining  species  is  adjusted  so  that  the  species  gradients  sum  to  zero  as  shown 

T  species  — 


N3 

\  ^ 


(VTfc)  -  VT, 


max  species 


(2.29) 


lk=l 


An  additional  issue  that  can  arise  is  when  one  LES  cell  has  a  small  non-zero  species  mass  value  and 
the  neighboring  cell  has  a  larger  value  that  will  cause  the  profile  to  go  below  zero  when  the  profile 
is  linearly  extrapolated.  In  this  case  the  gradient  of  the  species  with  the  small  mass  fraction  could 
again  be  set  to  zero  but  that  would  unnecessarily  eliminate  the  gradient  completely.  Instead,  all 
of  the  gradients  are  equally  scaled  so  that  any  species  that  contains  negative  mass  fractions  within 
the  profile  is  adjusted  to  have  a  minimum  value  of  zero  at  the  furthest  node  along  the  line.  All  of 
the  species  gradient  are  multiplied  by  a  scaling  factor  as  shown 

2  Y 

^ 1  mtn 


VY  =  VT 


VIA 


(2.30) 


This  method  meets  the  requirements  of  equations  2.27  and  2.28  while  maintaining  as  much  of  the 
original  gradients  as  possible.  A  similar  scaling  process  is  applied  to  the  temperature. 


To  illustrate  this  process,  the  initial  species  mass  fraction  profiles  for  an  arbitrary  four  species 
fluid  are  considered.  Figure  2.4  shows  two  cells  with  their  corresponding  species  mass  fractions  and 
temperature  values  that  were  calculated  at  the  resolved  scale.  The  first  step  of  the  initialization 
process  of  the  left  cell  aligns  the  LEM  line  with  the  maximum  temperature  gradient  and  calculates 

the  species  mass  fraction  gradients  along  this  line.  These  gradients  are  used  to  linearly  extrapolate 
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the  species  profiles  for  the  sub-grid  LEM  and  are  shown  in  Figure  2.5a.  From  the  plot  it  can  be 
seen  that  both  Yj  and  I4  have  negative  species  mass  fraction  values  for  a  portion  of  their  profiles. 
Since  I4  has  a  value  of  zero  at  the  cell  center,  as  shown  in  Figure  2.4,  the  gradient  must  be  zero 
to  avoid  a  negative  mass  fraction.  Additionally,  the  gradient  for  Y\  is  modified  so  that  the  sum  of 
the  gradients  add  up  to  zero  per  equation  2.28.  Species  I3  has  a  non-zero  value  at  the  cell  center 
but  the  gradient  is  large  enough  such  that  a  negative  species  mass  fraction  will  exist  at  the  extent 
of  the  domain.  This  issue  is  mitigated  by  scaling  all  the  profiles  so  that  I3  is  zero  at  the  far  left 
node.  Figure  2.5b  shows  the  fiual  initial  profile  for  this  LEM  line. 


Yi=  0.55 

Yi=  0.30 

Y2=  0.40 

Yz=  0.20 

Yj=  0.05 

Yj=  0.40 

Y<=  0.00  . 

Yt=  0.10  . 

T  =  300  K 

T  =  1500  K 

Figure  2.4: 


Adjacent  LES  cells  with  their  corresponding  mass  fractions  and  temperature. 
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(a)  Initial  species  mass  fraction  along  a  LEM  line  (b)  Re-scaled  species  mass  fraction  along  a  LEM  line 
with  negative  mass  on  the  left  side.  without  negative  mass. 


Figure  2.5:  Initialization  process  along  LEM  line. 


2.4.3  Numerical  Implementation  of  LEM  Processes 


In  practice  equations  2.18  and  2.19  are  solved  in  three  explicit  steps,  with  each  step  corresponding 
to  the  diffusion,  production,  and  stirring.  The  LEM  equations  are, 


dYk 

dt 


Hr  T 


ld_ 


P 


(Dk„ 


,m  d 


( WmixYk )  -  VfYk  + 


uk 


(2.31) 


dT 

~dt 


Ft,  stir  T 


pYpr 


ds  V  ds 


dT 

+  pTs 


„  Dk)M  d 

Cp^k  „  (IF, 

wmix  ds 


rYk) 


+ 


L jJt 


pCp.r 


(2.32) 


The  species  production  equations  are, 


dY^ 

dt 

dT 

~dt 


P 

Cjt 

pCp^mix 


(2.33) 

(2.34) 


An  ODE  solver  is  implemented  to  advance  equations  2.33  and  2.34  in  time  because  the  chemical  time 
scales  can  be  very  small.  The  stirring  process  is  an  instantaneous  rearrangement  event  (detailed  in 
section  2.2)  and  therefore  no  ODE  is  necessary. 
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The  LEM  solution  is  advanced  to  the  LES  physical  time  step  by  a  sub-grid  time  step  interval  of 
At,  as  seen  in  Figure  2.6.  The  value  of  At  is  determined  as  the  minimum  of  the  stirring  or  diffusion 
time  steps  (ref.  sections  2.4.3. 1  and  2. 4. 3. 2).  In  practice  At  typically  corresponds  to  the  stirring 
time  step  because  it  is  smaller  than  the  diffusion  time  step. 

The  combustion  equations  2.33  and  2.34  are  simultaneously  integrated  over  an  interval  of  At 
using  an  ODE  solver.  Next,  the  updated  solution  is  instantaneous  stirred  using  the  triplet  map 
rearrangement  process.  Finally,  the  diffusion  is  advanced  by  At  using  a  finite  difference  scheme 
which  smooths  out  the  sharp  edges  from  the  stirring.  This  process  is  repeated  until  the  LEM 
solution  has  matched  the  LES  interval,  with  a  partial  step  implemented  at  the  end  as  needed. 


_ 1 _ 1 _ 1 _ 1 _ 1 _ 1 

1  1  1  1  1  1  1  1 

At  Fractional  time  step 

as  necessary 

A  4- 

~ LES  H 

Figure  2.6:  LEM  time  stepping  schematic. 


2. 4. 3.1  Stirring 


The  stirring  is  implemented  using  three  parameters:  the  local  eddy  size,  £,  the  frequency,  /,  and 
the  stirring  location,  so-  The  eddy  size  is  chosen  randomly  from  a  probability  density  function 
(PDF),  /(£),  in  the  range  of  A  to  rj.  The  eddy  size  PDF  is: 

-f— 8/3 

/«)  =  „-S/!  _  A-S/3  <2-35> 

This  PDF  is  integrated  to  obtain  a  cumulative  distribution  function  (CDF).  Enforcing  the  condition 
that  the  CDF  takes  a  value  of  unity  when  £  =  A  results  in  the  following, 


CDF  (  A"5/3  -  r/-5/3 )  +  t? 


-5/3 


(2.36) 


The  length  of  the  eddy  for  each  time  step  is  obtained  by  choosing  a  random  number  for  the  CDF 
value  in  equation  2.36.  The  stirring  frequency,  /,  is  calculated  as: 


54  uReA  [(A/??)5/3  -  1 
ITCaA3  1-(A/t/)4/3 


(2.37) 
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where  C\  is  a  constant  determined  to  be  0.067  [16].  This  value  is  then  used  to  calculate  the  stirring 
time  step  as  follows: 

At  stir  =  j ^  (2-38) 

Once  the  eddy  size  is  calculated  it  is  then  possible  to  determine  the  location  of  the  stirring  event. 
The  location  is  chosen  using  a  random  number  generator,  assuming  an  initial  distribution,  with  the 
requirement  that  the  whole  eddy  length,  £,  must  to  fit  on  within  the  domain.  This  is  accomplished 
by  putting  a  restriction  that  the  node  at  which  the  event  occurs  must  be  less  than  IVlem  —  £■  Once 
the  eddy  length  and  location  are  determined,  the  triplet  map  can  then  be  applied  to  the  scalar 
(species  mass  fraction  and  temperature)  profile.  The  triplet  map  rearranges  the  profile  as  follows: 

0(s)  ->  0(M(s))  (2.39) 


where  the  scalar  value  is  mapped  to  the  one- dimensional  coordinate,  s,  and  M(s )  is  the  mapping 
[24]  which  is  defined  by: 


M(s) 


3  (s  -  s0) 


so  +  s 


2  i  -  3(s  -  s0) 
3(s  -  s0)  -  21 


s  -  s0 


if  so  <  s  <  so  +  1/3  i 
if  s0  +  1/3*  <  s  <  s0  +  2/3  t 
if  so  +  2/3£  <  s  <  so  +  £ 
if  otherwise 


2. 4. 3. 2  Diffusion 


The  diffusion  equations  in  LEM  are  solved  using  forward  Euler  in  time  and  second-order  central 
differencing  in  space.  The  equations  are: 


dYk  _  Id 
dt  p  ds 


P 
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dT 

~dt 


d  fdT\ 
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hbjnia;  ds 
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(' WmixYk )  -  VfYk 
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ds 


\k= 1 


Dk,M  d 

Wrrtir  ds 


(WmixYk) 


(2.40) 

(2.41) 


The  derivatives  in  the  above  equations  are  evaluated  using  central  differencing  stencil  as  shown: 

d(f>  </>i+ 1  —  (f>i- 1 


ds  2As 

A  Neumann  boundary  condition  =  0^  is  implemented  as: 
dcj)  —3(/»i  +  402  -  03 


=  0  =>■  0!  = 


402  -  03 


ds  2As 

<90 30jvz,BM+i  —  4 4>Nlem  4>Nj 

ds  2As 


LEM- 1  _  n  .  ±  _  4 0JVj 

—  _  U  ^  <PNlem  +  1  - 


LEM  ^Nl 

3 


f— i 


(2.42) 

(2.43) 

(2.44) 
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The  time  step  for  diffusion  is  chosen  as 


Atdifr  =  VNN - 

max 


As2 


where  VNN  is  set  at  0.5  for  numerical  stability. 


(2.45) 


2. 4. 3. 3  Combustion 


The  change  in  species  mass  fraction  and  temperature  due  to  combustion  is  calculated  by  the 
following  ODE’s 


dYk 

at 

dT 

~dt 


Uk 

P 

UJj' 

pCp^mix 


(2.46) 

(2.47) 


The  production  rate  for  each  species  is  expressed  as: 

M 

uJk  =  wk  ("k,j  -  ul,j)  wj  (2-48) 

3= 1 

where  M  is  the  number  of  reactions,  Wk  is  the  molecular  weight,  v'k  ■  and  v'k-  are  the  stoichiometric 
coefficients  of  the  reactants  and  products,  respectively,  and  Wj  is  the  progress  rate  variable.  The 
progess  rate  variable  can  be  calculated  for  global  reactions,  equilibrium  reactions,  Lindemann  falloff 
reactions,  third  body  reactions,  troe  reactions  and  Tsang-Herron  reactions.  The  ODE’s  are  solved 
using  the  solver,  DVODE,  which  uses  implicit  Adams  method  to  solve  equations  (2.46)  and  (2.47). 


2.4.4  Filtered  Species  Production  Term 


The  LEM  solution  advances  the  stirring,  diffusion,  and  species  production  to  the  next  physical 
time  step,  A f°r  the  entire  LEM  line.  The  filtered  species  production  term  for  each  species  is 
calculated  by  determining  the  average  species  production  over  the  entire  LEM  line  as  follows, 


1 

A'lem  +  1 


M.em+1 

i— 1 


(2.49) 


where  k  represents  each  individual  species. 
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2.5  Comparison  to  Past  Work 


This  approach  uses  LEM  to  obtain  the  filtered  species  production  term  in  the  large  scale  species 
mass  fraction  equations  2.4.  Previous  work  did  not  solve  the  Eulerian  form  of  the  species  equations 
at  the  resolved  scale  [43].  Instead,  the  species  mass  fraction  is  accounted  for  by  a  Lagrangian  splicing 
method  in  which  mass  is  transferred  between  LEM  lines  using  a  large  scale  advection  approach  as 
discussed  in  section  1.1.3.  The  large  scale  advection  is  an  ad  hoc  process  that  arbitrarily  transfers 
mass  between  LES  cells  accounting  only  for  advection,  while  neglecting  diffusion.  The  benefit  of 
the  current  approach  is  that  the  LEM  solution  contributes  to  the  full  multi-dimensional  species 
equations  at  the  resolved  scale,  which  includes  a  consistent  form  of  the  resolved  scale  convection 
and  diffusion. 

2.6  Assumptions  and  Limitations  of  LEM 

At  this  point  it  is  important  to  explicitly  state  the  assumptions  that  were  used  in  deriving  the 
species  transport  and  temperature  equations  to  better  understand  the  limitations  of  LEM.  A  salient 
feature  of  the  LEM  approach  is  that  a  one-dimensional  line  is  embedded  within  a  multi-dimensional 
domain.  It  is  aligned  with  the  maximum  temperature  gradient  in  order  to  capture  the  flame  front, 
however,  it  is  not  able  to  capture  all  of  the  relevant  flow  properties  that  would  be  captured  by  a 
three-dimensional  DNS.  Second,  the  form  of  the  energy  equation  used  in  LEM  assumes  that  there 
is  a  constant  pressure  within  the  LES  cell  from  one  physical  time  step  to  the  next.  This  makes 
LEM  appropriate  for  deflagarations  but  is  not  correct  in  areas  where  compressibility  effects  are 
important,  such  as  detonations  or  acoustics  coupling.  Third,  there  are  two  separate  temperatures 
calculated  because  an  energy  equation  is  solved  at  both  the  resolved  and  unresolved  scales.  The 
temperature  obtained  from  LEM  is  discarded  at  the  end  of  each  physical  time  step  since  it  is  more 
approximate.  Additionally,  a  no  flux  boundary  condition  is  applied  to  the  LEM  line  so  there  is  no 
transfer  of  mass  between  LEM  lines  through  diffusion  or  advection.  Lastly,  in  areas  where  the  LES 
grid  spacing  approaches  the  Kolmogorov  scales,  the  LEM  grid  will  reduce  to  a  single  node  or  a  few 
nodes.  In  this  limit,  this  LEM  approach  will  naturally  reduce  to  a  laminar  combustion  source  term 
closure  which  is  the  correct  DNS  limit. 
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CHAPTER  3 


Stand-Alone  LEM  Implementation 

The  LEM  equations  are  solved  in  three  explicit  processes,  which  includes  stirring,  diffusion  and 
combustion.  It  is  crucial  that  all  of  these  process  are  conservative  so  that  the  LEM  model  is 
physically  accurate.  This  chapter  tests  the  predictions  of  the  LEM  model  for  a  single  sub-grid 
line  and  then  focuses  on  verifying  that  the  model  for  the  three  processes.  In  addition,  the  energy 
spectrum  is  calculated  in  order  to  verify  that  LEM  can  accurately  represent  turbulent  statistics  in 
the  inertial  sub-range. 

3.1  Stirring 

The  stand-alone  sub-grid  stirring  process,  as  described  in  section  3.3.3. 1,  is  tested  using  four  species 
distributed  along  a  36  cell  LEM  line.  The  cell  centered  species  mass  fraction  values  for  Y\ ,  >2,  I3, 
and  I4  are  0.4,  0.3,  0.2,  and  0.1,  respectively,  which  could  represent  mass  fractions  of  oxidizer, 
products,  and  fuel,  for  example.  An  equal  and  opposite  gradient  is  applied  to  species  Y\  and  Y2, 
similar  to  Y3  and  Y4,  so  that  the  total  species  mass  fraction  within  the  domain  adds  to  unity  at 
each  LEM  node,  as  seen  in  Figure  3.1a.  Additionally,  a  cell  centered  temperature  value  of  300  K 
is  extrapolated  to  created  a  temperature  gradient  within  the  domain,  and  is  shown  in  Figure  3.1b. 
A  single  stirring  event  is  applied  at  a  random  location  along  the  LEM  line  (ref:  Figures  3.2a  and 
3.2b)  to  simulate  the  effect  of  sub-grid  advection  within  the  domain.  It  is  worth  noting  that  the 
same  eddy  size  and  location  is  applied  to  both  the  species  mass  fraction  and  temperature  profiles 
so  as  to  maintain  the  physical  representation  of  an  eddy  acting  on  the  sub-grid,  as  described  in 
section  2.2.  The  single  stirring  event  is  able  to  conserve  the  scalar  quantities  contained  in  the  LEM 
line  but  is  important  that  this  process  remains  conservative  for  a  large  number  of  instances. 

To  further  test  the  stirring  process,  the  code  is  run  for  one  hundred  stirring  events  (ref:  Figures 
3.3a  and  3.3b  )  on  the  same  initial  profiles  and  the  total  mass  of  each  species  and  temperature 


(a)  Initial  species  mass  fraction  profile.  (b)  Initial  temperature  profile. 

Figure  3.1:  Initial  species  mass  fraction  and  temperature  distribution  for  a  single  LEM  line. 

contained  within  the  LEM  line  is  calculated  before  and  after  the  stirring  using  a  trapezoidal  rule, 
as  shown  in  Table  3.1.  The  results  in  the  table  show  that  both  the  species  mass  fraction  and 
temperature  quantities  are  conserved  for  a  large  number  of  stirring  events,  which  verifies  that  the 
stirring  process  is  correctly  implemented.  In  practice,  a  typical  LES  time  step  will  only  have  a 
handful  of  sub-grid  stirring  events  which  means  that  this  test  far  exceeds  the  number  of  stirring 
events  that  would  be  applied  during  a  given  simulation. 
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(a)  Single  stirring  event  applied  to  the  species 
mass  fraction  profiles. 


file. 


Figure  3.2:  Single  stirring  event  applied  to  the  species  mass  fraction  and  temperature  profiles. 


(a)  One  hundred  stirring  event  applied  to  the  species 
mass  fraction  profiles. 


(b)  One  hundred  stirring  event  applied  to  temperature 
profile. 


Figure  3.3:  One-hundred  stirring  events  applied  to  the  species  mass  fraction  and  temperature 
profiles. 
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Table  3.1:  Area  before  and  after  100  stirring  events 


Variable 

Area  Before  Stirring 

Area  After  Stirring 

V 

4.00  xlO”4 

4.00  xl0~4 

3.00  xlO-4 

3.00  xlO”4 

^3 

2.00  xl0~4 

2.00  xlO-4 

n 

1.00  xlO-4 

1.00  xlO-4 

total  mass 

1.00  xlO"3 

1.00  xlO-3 

T(K) 

300.0 

300.0 

3.2  Diffusion 

The  diffusion  process,  as  described  in  Section  2. 4. 3. 2,  is  tested  using  four  species  distributed  along  a 
36  cell  LEM  line.  The  same  initial  distributions  as  the  previous  section  are  used  here  and  are  again 
shown  in  Figures  3.4a  and  3.4b.  The  diffusion  process  is  run  with  a  Neumann  boundary  condition 
for  a  total  of  200  ms  which  corresponds  to  around  200  K  time  steps.  The  no-flux  boundary  condition 
is  implemented  so  the  total  mass  and  temperature  contained  within  the  LEM  line  is  conserved. 

The  diffused  profiles  are  shown  in  Figures  3.5a  and  3.5b  and  it  can  be  seen  that  all  the  profiles 
have  dissipated  close  to  the  nominal  values,  as  expected.  The  quantity  of  species  and  the  tempera¬ 
ture  contained  in  the  domain  is  calculated  before  and  after  the  diffusion  process  using  a  trapezoidal 
rule  and  are  shown  in  Table  3.2.  The  tabulated  results  show  that  the  diffusion  process  has  ac¬ 
cumulated  error  that  varies  between  1%  and  3%  of  the  expected  values.  This  error  is  due  to  the 
second  order  finite  difference  scheme  that  is  used  to  represent  the  spatial  gradient  and  eliminating 
this  error  completely  is  unavoidable.  However,  in  practice  the  diffusion  will  take  place  for  a  much 
shorter  amount  of  time  and  for  fewer  time  steps  than  in  this  test,  so  the  accumulated  error  with 
this  process  will  be  minimized  in  a  given  simulation. 
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Figure  3.4:  Initial  species  mass  fraction  and  temperature  distribution  for  a  single  LEM  line. 
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(a)  Diffused  species  mass  fraction  profiles. 


(b)  Diffused  temperature  profile. 


Figure  3.5:  Diffusion  process  applied  to  the  species  mass  fraction  and  temperature  profiles. 
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Tabic  3.2:  Area  before  and  after  diffusion. 


Variable 

Area  Before  Diffusion 

Area  After  Diffusion 

V 

4.00  xlO"4 

4.03  xl0~4 

>"2 

3.00  xl0“4 

2.96  xl0~4 

y3 

2.00  xl0~4 

2.04  xlO-4 

y4 

1.00  xHT4 

0.97  xHT4 

total  mass 

1.00  xlO-3 

1.00  xlO-3 

T(K) 

300.0 

291.12 

3.3  Combustion 

The  stand-alone  combustion  process  was  tested  using  GRI-Mech  1.2  [44],  GRI-Mech  is  an  optimized 
detailed  chemical  reaction  mechanism  capable  of  the  best  representation  of  natural  gas  flames  and 
ignition  to  date.  GRI-Mech  is  a  list  of  elementary  chemical  reactions  and  associated  rate  constant 
expressions.  The  version  utilized  in  this  work  contains  177  reactions  and  32  species  which  is 
commonly  used  because  of  its  relatively  small  size  and  ability  to  accurately  predict  methane  and 
air  based  combustion. 

A  simulation  is  run  using  LEM  combustion  and  it  is  compared  to  the  solution  from  CHEMKIN, 
a  software  tool  developed  by  Sandia  National  Laboratories  for  solving  complex  chemical  kinetics 
problems;  the  mechanism  in  GRI-Mech  1.2  was  incorporated  into  CHEMKIN  as  well.  The  simula¬ 
tion  is  initialized  with  species  mass  fraction  values  of  0.504  for  oxygen  (O2),  0.370  for  water  (H2O), 
and  0.126  for  methane  (CH4)  at  101  kPa  and  1500  K.  It  can  be  seen  from  Figure  3.6  that  both 
methods  achieve  the  same  species  mass  fraction  values  for  oxygen  and  methane  but  in  a  slightly 
different  amount  of  time.  The  combustion  temperature  is  compared  in  Figure  3.7  and  again  shows 
very  good  agreement  in  final  temperature  but  has  the  same  discrepancy  in  the  reaction  time.  The 
difference  in  time  is  probably  due  to  the  lower  tolerance  value  set  in  CHEMKIN  compared  to  the 
LEM  solver. 
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Figure  3.6:  Comparison  of  oxygen  and  methane  during  combustion  for  GRI-Mech  1.2  reaction 
mechanism  using  LEM  and  CHEMKIN. 


Figure  3.7:  Comparison  of  temperature  during  combustion  for  GRI-Mech  1.2  raction  mech¬ 
anism  using  LEM  and  CHEMKIN. 
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3.4  Energy  Spectrum  of  LEM 


The  kinetic  energy  contained  in  a  turbulent  flow  cascades  down  from  the  largest  to  the  smallest 
eddies.  Since  turbulence  contains  a  continuous  spectrum  of  scales,  it  is  convenient  to  do  analysis  in 
terms  of  the  spectral  distribution  of  energy  [5].  The  wavenumber  spectrum  in  turbulent  flows  can 
be  divided  into  three  sections,  which  are  the  energy  containing  eddies,  the  inertial  sub-range,  and 
the  viscous  range.  In  LES,  the  resolved  scale  captures  the  energy  containing  eddies  and  extends 
into  the  inertial  range.  The  LEM  sub-grid  model  aims  to  represent  the  inertial  sub-range  of  a 
turbulent  flow  which  extends  down  to  the  Kolmogorov  length  scale. 

A  well  established  property  of  the  inertial  sub-range  is  the  Kolmogorov  -5/3  law.  It  is  important 
that  the  LEM  be  compared  to  this  criterion  in  order  to  validate  its  viability  as  a  turbulence  model. 

The  energy  spectrum  produced  by  the  LEM  model  is  determined  by  running  the  full  LEM  solver, 
as  described  in  section  2.4.3,  using  a  single-step  methane  reaction  that  contains  four  different  species 
within  the  domain.  The  stoichiometric  reaction  is 

CH4  +  202  =  2H20  +  C02  (3.1) 

The  two  cases  that  are  considered  here  are  a  5  mm  long  LEM  line  that  contains  96  LEM  cells  and  2 
mm  long  LEM  line  that  contains  120  LEM  cells.  These  two  cases  are  chosen  because  they  represent 
possible  LEM  grids  in  practical  simulations.  Case  1  (5mm)  has  and  initial  cell  centered  species  mass 
fraction  values  of  0.4,  0.3,  0.2,  and  0.1  for  oxygen,  methane,  water,  and  carbon  dioxide,  respectively. 
Case  2  (2mm)  has  the  stoichiometric  values  of  0.8  for  oxygen  and  0.2  for  methane,  hence  without 
the  presence  of  product  at  time  t=0.  Both  cases  are  run  for  a  total  of  10  fj, s,  and  the  resulting 
stirred  profiles  for  methane  are  shown  in  Figures  3.8  and  3.10,  for  case  1  and  2,  respectively. 

To  examine  the  frequency  content  of  the  signal,  a  power  spectral  density  analysis  is  performed 
on  the  final  profiles.  A  total  of  5  mm  of  spatial  data  are  analyzed,  with  the  data  spaced  52  /.ini 
apart  for  case  1.  For  case  2  a  total  of  2  mm  of  spatial  data  are  analyzed,  with  the  data  spaced  16 
gm  apart.  The  PSD  results  are  plotted  on  a  log  scale  and  are  shown  in  Figures  3.9  and  3.11,  for 
cases  1  and  2,  respectively.  Both  figures  have  an  additional  red  line  with  a  -5/3  slope  that  is  used 
verify  that  the  LEM  turbulence  does  indeed  represent  the  inertial  sub-range. 
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Figure  3.8:  The  methane  profile  along  the  5  mm  long  domain  after  LEM  has  been  run  for 
10  /xs,  with  initial  conditions  corresponding  to  case  1. 


Figure  3.9:  PSD  of  the  methane  profile  for  the  5  mm  long  line  after  LEM  has  been  run  for 
10  /us,  with  initial  conditions  corresponding  to  case  1.  The  red  line  has  a  -5/3  slope  that  is 
used  to  verify  the  turbulence  profile  follows  Kolomogorov  -5/3  law. 
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Figure  3.10:  The  methane  profile  along  the  2  mm  long  domain  after  LEM  has  been  run  for 
10  n s,  with  initial  conditions  corresponding  to  case  2. 


Figure  3.11:  PSD  of  the  methane  profile  for  the  2  mm  long  line  after  LEM  has  been  run  for 
10  /us,  with  initial  conditions  corresponding  to  case  2.  The  red  line  has  a  -5/3  slope  that  is 
used  to  verify  the  turbulence  profile  follows  Kolomogorov  -5/3  law. 
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3.5  Summary  of  Stand-Alone  Results 


The  three  fundamental  processes  of  stirring,  diffusion,  and  chemical  reaction  associated  with  LEM 
were  tested  to  verify  that  each  performs  correctly.  The  stirring  successfully  implements  a  triplet  map 
to  produce  a  turbulent  profile  that  is  completely  conservative.  The  diffusion  functions  as  expected 
but  produces  a  small  error  due  to  truncation  errors  on  a  finite  mesh.  The  LEM  combustion  model 
was  tested  using  the  Gri-Mech  reaction  mechanism  and  produces  very  similar  results  to  Chemkin, 
noting  that  the  slight  difference  in  reaction  time  is  probably  a  result  of  the  time  stepping  method 
employed  in  the  Chemkin  solver. 

Additionally,  a  power  spectral  density  analysis  has  been  used  to  determine  the  energy  spectrum 
produced  by  the  LEM  model,  and  the  results  show  that  applying  LEM  to  a  scalar  profile  produces 
a  -5/3  slope  in  the  energy  spectrum  that  is  expected  to  be  seen  in  the  inertial  sub-range.  These 
results  validate  LEM  as  a  reasonable  turbulence  model  for  combustion  processes.  This  chapter 
has  verified  that  the  LEM  model  has  been  properly  implemented  and  appropriately  represents 
turbulence  in  the  sub-grid  scale.  The  following  chapter  will  couple  LEM  with  GEMS  to  form  the 
LEM-LES  model.  The  combined  solver  will  then  be  applied  to  a  bluff-body  stabilized  flame. 
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CHAPTER  4 


LEM-LES 

This  chapter  uses  the  integrated  LEM-LES  solver  to  simulate  reactive  flow  in  the  Volvo  rig,  which 
consists  of  a  rectangular  duct  that  contains  an  equilateral  triangle-shaped  bluff-body  flameholder. 
The  results  for  temperature,  centerline  axial  velocity,  normalized  velocity,  and  lean  blow-out  equiv¬ 
alence  ratio  from  the  reactive  LEM-LES  code  are  compared  to  experimental  data  to  assess  the 
validity  of  LEM  as  a  sub-grid  model  for  turbulent  combustion  in  bluff-body  stabilized  flames. 

4.1  Volvo  Rig  Description 

The  Volvo  rig  was  designed  in  the  early  1990s  and  has  been  a  valuable  experiment  for  a  compu¬ 
tational  comparison  due  to  the  detailed  experimental  data  available  [27,  28,  29].  The  test  rig  is 
a  1.55m  long  rectangular  duct  that  is  0.24m  wide  by  0.12m  high  and  is  fed  by  a  choked  air  inlet. 
The  rig  is  divided  into  a  0.55m  inlet  section  and  a  l.OOrn  combustor  section.  The  inlet  section 
injects  a  choked  air  and  fuel  mixture  into  the  duct  which  then  travels  through  a  honeycomb  and 
screen  before  it  enters  the  combustor  section.  The  combustor  section  contains  a  40mm  equilateral 
triangular  bluff  body  that  sits  0.318m  downstream  from  the  inlet  section  which  spans  the  width  of 
the  0.24m  duct  and  acts  as  a  flame  holder.  After  exiting  the  combustor  the  flow  is  exhausted  into 
a  large  circular  duct  that  is  0.352m  in  diameter.  Figure  4.1  shows  the  geometric  schematic  for  the 
Volvo  test  rig  downstream  of  the  honeycomb  screen. 

4.2  Simulation  Description 

The  LEM-LES  solver  is  applied  to  the  Volvo  rig  geometry  to  compute  lean  blow-out,  time  averaged 
temperature,  and  normalized  time  averaged  velocity.  The  simulation  is  performed  in  two  dimensions 
using  a  mesh  that  contains  49,978  cells.  The  downstream  region  of  the  mesh  is  important  because 
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Figure  4.1:  Cross-sectional  view  of  the  Volvo  geometry,  including  the  triangular  flameholder. 

that  is  where  the  flame  resides  and  therefore  it  is  important  to  properly  capture  this  region.  A  cell 
width  of  4nmi  is  chosen  in  the  axial  direction  downstream  of  the  flameholder  because  it  matches 
Pope’s  criterion  in  which  80%  of  the  turbulent  kinetic  energy  is  resolved. 

The  inlet  boundary  has  a  mass  inflow  of  0.62  kg/s  and  a  stagnation  temperature  of  288K.  The 
flameholder  and  channel  walls  are  treated  as  adiabatic,  non-slip  surfaces,  meanwhile  the  exit  has 
a  fixed  back  pressure  of  lOOkPa.  The  combustion  mechanism  is  a  global  propane  reaction  that 
contains  6  species:  propane  (C3H8),  oxygen  (O2),  nitrogen  (N2),  carbon  dioxide  (CO2),  water 
(H2O),  and  neon  (Ne).  The  simulations  are  run  with  a  physical  time  step  of  5xl0~5s  using  both 
LEM-LES  and  a  laminar  closure  for  the  combustion  source  term  so  that  the  results  can  be  compared. 
The  temperature  and  velocity  profiles  are  computed  for  an  equivalence  ratio  of  0.65  to  compare 
with  experimental  data  [28].  To  predict  lean  blow-out,  an  equivalence  ratio  of  0.5  is  chosen  because 
that  has  been  shown  in  experiments  to  extinguish  the  flame  [27]. 
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4.3  Results 


4.3.1  Stable  Combustion 

The  first  simulation  that  is  run  compares  laminar  reactive  flow  code  results  to  the  LEM-LES  solver 
with  an  equivalence  ratio  of  0.65.  A  single  flow  through  time  for  the  input  parameters  discussed  in 
Section  4.2  is  0.087  seconds.  The  simulations  were  run  for  0.15  and  0.25  seconds  which  corresponds 
to  roughly  2  and  3  flow  through  times,  respectively.  The  instantaneous  temperature  contour  plots 
are  shown  in  Figures  4.2  and  4.3  for  results  from  the  two  different  methods  at  times  0.15  and  0.25 
seconds,  respectively.  The  contour  plots  show  that  the  laminar  reactive  flow  solution  achieves  a 
temperature  around  1950K  downstream  of  the  flameholder  where  LEM-LES  predicts  1850K  for 
the  same  region.  A  qualitative  look  at  the  instantaneous  profiles  shows  that  the  contours  have 
noticeably  different  flame  shapes  downstream  of  the  flameholder  in  addition  to  the  quantitative 
temperature  difference.  The  instantaneous  contour  profiles  can  be  misleading  however,  because 
the  models  are  initialized  in  a  non-physical  manner  which  may  result  in  the  simulations  being  at 
different  points  in  the  combustion  cycle.  To  ameliorate  this  issue  an  average  temperature  profile 
is  useful  because  it  will  capture  the  whole  combustion  cycle  throughout  the  simulation.  The  time 
averaged  temperature  shown  in  Figure  4.4  reveals  that  the  downstream  profile  flame  shapes  are 
in-fact  similar  with  the  exception  of  the  temperature  value  downstream  of  the  flameholder,  with 
higher  temperatures  for  the  laminar  closure  model  persisting  into  the  farfield. 

A  quantitative  analysis  is  performed  by  time  averaging  the  temperature  and  normalizing  velocity 
profiles  downstream  of  the  flameholder  and  comparing  these  values  to  experimental  data.  The 
locations  that  are  measured  for  the  time  averaged  temperature  profile  are  0.15m,  0.35m,  and  0.55m 
downstream  of  the  flameholder.  The  normalized  time  averaged  velocity  is  measured  0.015m  and 
0.15m  downstream  of  the  flameholder.  In  addition  to  the  transverse  velocity  a  normalized  centerline 
axial  velocity  is  measured  starting  from  the  edge  of  the  flameholder  and  extending  the  length  of 
the  duct.  The  location  of  the  measured  data  is  shown  in  Figure  4.5. 

The  time  averaged  temperature  profiles  (Figures  4.6,  4.7,  and  4.8)  show  that  both  LEM-LES 
and  laminar  combustion  over-predict  the  temperature  in  all  three  locations.  The  CFD  solutions 
furthermore  become  more  erroneous  further  downstream  of  the  flameholder.  The  LEM-LES  based 
temperatures  are  lower  at  each  point  and  are  therefore  closer  to  the  experimental  data  than  tem- 
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Temperature,  K 
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(a)  Instantaneous  temperature  contours  of  laminar  closure  model  0.15  seconds  after  initialization 
with  an  equivalence  ratio  of  0.65. 


(b)  Instantaneous  temperature  contours  of  LEM  sub-grid  model  0.15  seconds  after  initialization 
with  an  equivalence  ratio  of  0.65. 

Figure  4.2:  Comparison  of  the  instantaneous  temperature  contours  of  laminar  combustion 
to  LEM  sub-grid  model  0.15  seconds  after  initialization  with  an  equivalence  ratio  of  0.65. 

peratures  derived  from  the  laminar  reactive  flow  calculation.  It  can  also  be  seen  that  both  of  the 
solutions  predict  the  width  of  the  flame  sufficiently  well  but  produce  a  steeper  gradient  at  the 
fringes  of  the  flame  than  is  seen  in  the  experimental  data. 
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Temperature,  K 
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(a)  Temperature  contours  of  laminar  closure  model  0.25  seconds  after  initialization  with  an  equiv¬ 
alence  ratio  of  0.65. 


(b)  Temperature  contours  of  LEM  sub-grid  model  0.25  seconds  after  initialization  with  an  equiv¬ 
alence  ratio  of  0.65. 

Figure  4.3:  Comparison  of  temperature  contours  of  laminar  combustion  to  LEM  sub-grid 
model  0.25  seconds  after  initialization  with  an  equivalence  ratio  of  0.65. 

The  time  averaged  normalized  velocity  profiles  (Figures  4.9  and  4.10)  show  that  the  LEM-LES 
and  laminar  combustion  predict  similar  velocity  profiles  at  both  locations.  The  LEM-LES  velocity 
profile  0.015m  downstream  of  the  flameholder  (Figures  4.9)  does  not  match  the  experimental  data 
along  the  centerline  as  well  as  the  laminar  reactive  flow  model,  however,  both  models  closely  match 
the  experimental  data  in  the  regions  away  from  the  centerline.  Figure  4.10  shows  that  0.15m 
downstream  of  the  flameholder  both  models  are  again  very  similar  and  over-predict  the  centerline 
velocity.  This  higher  centerline  velocity  prediction  occurs  for  both  models  in  accordance  with  the 
downstream  measured  locations,  as  shown  in  Figure  4.11.  It  can  be  seen,  as  in  the  temperature 
profiles,  that  the  LEM-LES  solver  overall  predicts  values  closer  to  the  experimental  data  than  does 
the  laminar  reactive  flow  model. 
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(a)  Time  averaged  temperature  contour  of  laminar  closure  model  with  an  equivalence  ratio  of  0.65. 


(b)  Time  averaged  temperature  contour  of  LEM  sub-grid  model  with  an  equivalence  ratio  of  0.65. 
Figure  4.4:  Comparison  of  time  averaged  temperature  contours  of  laminar  reactive  flow  (a) 
to  LEM  sub-grid  model  (b)  with  an  equivalence  ratio  of  0.65. 
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Figure  4.5:  Locations  downstream  of  the  flameholder  where  the  time  averaged  temperature 
and  velocity  prohles  are  measured. 
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Figure  4.6:  Temperature  profile  0.15m  downstream  of  flameholder  for  laminar  combustion 
(red  line)  and  LEM  sub-grid  (blue  line)  model  compared  to  experimental  data  (black  dots). 


Figure  4.7:  Temperature  profile  0.35m  downstream  of  flameholder  for  laminar  combustion 
(red  line)  and  LEM  sub-grid  (blue  line)  model  compared  to  experimental  data  (black  dots). 
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Figure  4.8:  Temperature  profile  0.55m  downstream  of  flameholder  for  laminar  combustion 
(red  line)  and  LEM  sub-grid  (blue  line)  model  compared  to  experimental  data  (black  dots). 


Figure  4.9:  Normalized  axial  velocity  profile  0.015m  downstream  of  flameholder  for  laminar 
combustion  (red  line)  and  LEM  sub-grid  (blue  line)  model  compared  to  experimental  data 
(black  dots). 


46 


1.5 


Q 


1.0  - 


0.5 


0.0  - 


-0.5 


-1.0  - 


-1.5. 


Laminar  combustion 
LEM  sub-grid  model 
Experimental  data 


-1.0  -0.5 


2.5 


Figure  4.10:  Normalized  axial  velocity  profile  0.15m  downstream  of  flameliolder  for  laminar 
combustion  (red  line)  and  LEM  sub-grid  (blue  line)  model  compared  to  experimental  data 
(black  dots). 
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Figure  4.11:  Centerline  axial  velocity  of  laminar  combustion  (red  line)  and  LEM  sub-grid 
(blue  line)  model  compared  to  experimental  data  (black  dots). 
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The  LEM  model  is  a  multi-scale  approach  that  solves  an  energy  equation  at  both  the  resolved 
and  sub- grid  scales  that  results  in  two  separate  temperatures  as  discussed  in  Section  2.6.  It  is 
important  to  verify  that  the  difference  between  these  two  temperatures  is  sufficiently  small  because 
an  assumption  of  LEM  is  that  the  temperature  can  be  discarded  at  the  end  of  each  physical  time 
step.  The  magnitude  of  the  difference  between  the  LEM  sub-grid  temperature  and  the  LES  resolved 
temperature  is  calculated  for  the  same  flow  parameters  described  in  Section  4.2  and  is  shown  in 
Figure  4.12.  The  results  show  that  there  is  only  a  difference  in  temperature  on  the  edges  of  the 
flame  and  the  maximum  is  4K. 


Temperature,  K 

0  0.5  1  1.5  2  2.5  3  3.5  4 


Figure  4.12:  Magnitude  of  temperature  difference  between  the  resolved  (LES)  and  sub-grid 
(LEM)  temperature  values. 

4.3.2  Lean  Blow-out 

In  order  to  test  if  the  LEM-LES  solver  will  properly  predict  lean  blow-out  the  simulation  was 
run  on  the  same  mesh  as  before  but  with  an  equivalence  ratio  of  0.50.  Experiments  have  shown 
that  the  flame  will  blow-out  at  an  equivalence  ratio  of  0.55  so  the  chosen  equivalence  ratio  should 
produce  blow-out.  The  simulation  was  initialized  by  running  an  equivalence  ratio  0.65  for  three 
flow  through  times  to  stabilize  the  flame.  The  inflow  condition  is  then  changed  to  an  equivalence 
ratio  of  0.50  and  the  simulation  is  continued.  Figure  4.13a  shows  the  instantaneous  temperature 
contour  of  the  stabilized  combustion  after  the  first  three  flow  through  times  with  an  equivalence 
ratio  of  0.65.  Figure  4.13b  shows  the  instantaneous  temperature  contours  five  flow  through  times 
after  the  equivalence  ratio  has  been  lowered  to  0.50.  It  can  be  seen  from  these  results  that  the 
flame  temperature  has  significantly  decreased  with  a  reduction  in  equivalence  ratio,  as  would  be 
consistent  with  leaner  combustion,  but  it  does  not  achieve  lean  blow-out.  The  temperature  directly 
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behind  the  flame  holder  has  decreased  from  1780  K  to  1525  K  and  has  re-stabilized  at  the  lower 


temperature  as  seen  in  Figure  4.14. 
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(a)  Instantaneous  temperature  contour  using  the  LEM-LES  solver  0.26  seconds  after  initialization 
with  an  equivalence  ratio  of  0.65. 


(b)  Instantaneous  temperature  contour  using  the  LEM-LES  solver  0.412  seconds  after 
changing  the  equivalence  ratio  from  0.65  to  0.50. 


Figure  4.13:  LEM-LES  solver  applied  to  the  Volvo  rig  to  predict  lean  blow-out. 


4.4  Discussion 

The  LEM  model  is  applied  to  the  Volvo  experimental  rig  to  assess  its  viability  as  a  sub-grid  model 
for  an  LES  solver.  The  results  show  that  LEM-LES  solver  reasonably  matches  the  experimental 
temperature  downstream  of  the  flameholder  with  an  equivalence  ratio  of  0.65.  The  LEM-LES 
temperature  is  closer  to  the  experimental  data  than  a  comparable  laminar  combustion  model. 
LEM  over-predicts  the  centerline  velocity  but  again  comes  closer  to  the  experimental  data  than 
the  laminar  solution.  The  addition  of  diffusion  and  advection  of  the  LEM  model  have  a  noticeable 
effect  on  the  solution.  The  difference  in  LES  and  LEM  temperature  is  determined  to  be  no  greater 
than  4  K  at  any  point  in  the  simulation.  This  result  means  that  there  is  no  significant  inconsistency 
in  the  resolved  scale  and  sub-grid  representation  of  the  energy  equation. 

The  LEM-LES  solver  does  not  properly  predict  lean  blow-out  at  the  experimentally-determined 
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Figure  4.14:  Temperature-time  history  after  the  equivalence  ratio  has  been  changed  from 
0.65  to  0.50  at  a  point  2.5  cm  downstream  of  flameholder  and  1.5  cm  off  the  center  line. 

equivalence  ratio.  The  flame  has  been  shown  in  experiments  to  blow-out  at  an  equivalence  ratio  of 
0.55.  In  this  work  the  LEM-LES  solver  does  not  capture  blow-out  even  with  a  lower  equivalence 
ratio  of  0.50.  The  results  show  that  the  flame  temperature  behind  the  flameholder  has  decreased 
by  255K  but  has  reached  a  relatively  constant  temperature,  continuing  to  sustain  combustion. 
The  flame  however  is  starting  to  show  signs  of  blowing  out  in  the  region  far  downstream  of  the 
flameholder.  This  can  be  seen  in  the  Karrnan  vortex  shedding  type  profile  of  temperature  that  is 
present  in  this  region  (Figure  4.13b).  It  had  been  expected  that  the  LEM  model  would  predict  lean 
blow-out  at  an  equivalence  ratio  closer  to  the  experimental  data  because  the  flame  temperature  is 
lower  for  a  given  equivalence  ratio  than  with  the  comparable  laminar  combustion  closure  model. 

A  possible  reason  for  the  discrepancy  is  with  the  initialization  of  the  LEM  profiles  at  the 
beginning  of  each  sub-grid  calculation  (ref  Section  2.4.2).  At  each  LES  time  step  the  LEM  solution 
is  re-initialized  to  a  linear  profile.  This  means  that  the  turbulence  from  earlier  sub-grid  calculations 
(i.e.  previous  LES  time  step)  is  neglected.  The  model  only  performs  stirring  for  a  small  number 
of  times  during  each  sub-grid  calculation  which  could  mean  that  the  model  does  not  produce  the 
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appropriate  degree  of  turbulence.  The  lack  of  proper  representation  of  turbulence  is  one  possible 
reason  for  the  improper  prediction  of  lean  blow-out. 

Additionally,  the  reaction  mechanism  used  in  this  simulation  is  a  single  step  propane  reaction. 
This  is  a  significant  simplification  when  compared  to  a  physical  flame  that  contains  hundreds 
of  associated  reactions  in  conjunction  with  the  global  reaction.  This  current  approach  does  not 
represent  processes  such  as  the  formation  of  hydrogen  radicals  that  are  generally  involved  with 
ignition  and  extinction  and  can  cause  the  lean  blow-out  prediction  to  be  erroneous. 
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CHAPTER  5 


Effects  of  Boundary  Treatment  on  Thermo-acoustic 

Predictions 

Numerical  simulations  are  is  an  invaluable  tool  for  engineers  because  they  can  be  used  to  optimize 
and  eliminate  design  configurations  before  test  articles  are  built.  Simulations  can  also  provide  in¬ 
sight  and  diagnostics  for  regions  which  are  difficult  to  experimentally  interrogate.  However,  CFD 
calculations  can  take  months  of  computational  time  to  produce  a  solution  for  complex  unsteady 
configurations.  For  bluff-body  stabilized  flames,  a  common  simplification  used  to  reduce  com¬ 
putational  time  is  to  eliminate  the  farfield  portion  of  the  plenum  or  exhaust  section  behind  the 
combustor.  The  risk  in  reducing  the  size  of  the  domain  is  that  the  artificial  boundary  conditions 
that  are  imposed  may  not  be  physically  representative  of  the  actual  experiment,  and  as  the  bound¬ 
ary  moves  closer  to  the  region  of  interest  it  may  alter  the  solution.  This  chapter  focuses  on  the 
effect  of  the  outlet  boundary  condition  on  the  thermo-acoustic  response  inside  a  combustor.  The 
experiment  chosen  for  this  study  is  a  bluff  body  rig  used  at  the  AFR.L.  CFD  is  used  to  examine 
the  effect  of  multiple  outlet  boundary  conditions  in  an  effort  to  produce  an  optimum  design  for  the 
experiment  while  simultaneously  showing  the  effect  that  the  outlet  boundary  condition  can  have 
on  the  modeled  flow  inside  the  combustor. 

5.1  Description  of  the  Bluff-Body  Configuration 

The  rig  used  in  this  study  was  designed  by  AFRL  at  Wright-Patterson  Air  Force  Base  to  study  bluff 
body  stabilized  turbulent  flames  to  better  understand  the  behavior  of  augmentors  in  gas  turbine 
engines.  The  test  section  is  a  12.7  cm  (5  in)  by  15.24  cm  (6  in)  rectangular  duct  with  a  3.81  cm  (1.5 
in)  bluff  body  equilateral  triangle  spanning  the  15.24  cm  width.  Premixed  fuel  and  air  is  injected 
through  a  perforated  plate  upstream  of  the  test  section  with  flow  rates  ranging  from  0.181  kg/s  (0.4 
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lbm/s)  to  0.318  kg/s  (0.7  lbm/s).  The  flow  exits  the  duct  and  travels  into  a  34.70  cm  (13.66  in) 
circular  exhaust  tube  60.96  cm  (24  in)  downstream  of  the  duct  exit.  Figure  5.1  shows  the  schematic 
of  the  rig. 


(39.6  in)  (22.08  in)  (24.0  in)  (15.50  in) 

Figure  5.1:  Cross-sectional  view  of  the  bluff-body  combustor  rig  geometry. 

Within  practical  flow  regimes,  the  triangular  flameholder  in  the  rig  creates  a  recirculation  zone 
directly  downstream  of  its  trailing  edge,  trapping  hot  products  and  serving  as  a  source  of  heat  to 
ignite  the  incoming  fuel  and  oxygen  mixture.  A  shear  layer  forms  downstream  from  the  flameholder 
between  the  main  flow  and  recirculation  region,  allowing  for  improved  mixing  between  the  hot  and 
cold  gases.  Figure  5.2  shows  the  typical  re-circulation  region  created  by  the  flameholder  for  the 
flow  conditions  studied  here. 


Figure  5.2:  Recirculation  behind  a  triangular  bluff  body  (Adapted  from  [45]). 
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5.2  Preliminary  Design  Study 


A  concern  for  the  experimental  rig  designers  was  that  air  from  the  downstream  farfield  would  get 
drawn  into  the  test  section  and  influence  the  flowfield.  The  addition  of  a  nozzle  to  the  downstream 
end  of  the  rig  was  proposed  to  minimize  this  effect.  The  higher  speed  exhaust  would  prevent  any 
mass  from  the  surroundings  from  being  drawn  into  the  test  section.  The  idea  was  explored  by 
running  a  computational  study  that  added  a  nozzle  to  the  duct  exit.  The  streamlines  and  the 
pressure  contours  were  examined  to  see  how  the  outside  air  interacted  with  the  exhaust  air  from 
the  duct  exit.  The  study  was  performed  with  three  configurations;  a  baseline  configuration  with 
no  nozzle,  a  5°  nozzle,  and  a  10°  nozzle  at  the  duct  exit.  These  configurations  are  shown  in  Figure 
5.3.  The  CFD  simulations  were  conducted  in  GEMS  using  a  two-dimensional  non-reacting  fluid 
flow  to  study  the  effect  of  the  nozzle.  The  flow  was  driven  by  injecting  hot  air  at  2200  K  into  the 
duct.  The  farfield  was  initialized  with  atmospheric  conditions,  air  at  a  pressure  of  101  kPa  and  a 
temperature  of  288  K.  The  simulations  were  run  for  three  flow-through  times  so  that  the  results 
could  be  analyzed  after  the  initial  transients  had  dissipated. 

A  comparison  of  the  vorticity  contours  resulting  from  these  configurations  shows  that  all  three 
cases  exhibit  Karman  vortex  shedding  off  of  the  flameholder,  as  seen  in  the  right  side  of  Figure  5.4. 
The  presence  of  the  nozzle  does  not  appear  to  affect  this  phenomenon  in  any  appreciable  way.  The 
pressure  in  the  duct,  however,  is  sensitive  to  the  presence  of  the  nozzle,  as  indicated  in  the  left  side 
of  Figure  5.4.  It  can  be  seen  that  as  the  nozzle  angle  increases,  so  does  the  bulk  pressure  in  the 
duct.  This  can  easily  be  seen  in  Figure  5.4c  where  the  scale  is  adjusted  to  capture  the  pressure 
variation  in  the  vortex  shedding.  In  this  configuration  the  maximum  pressure  is  104  kPa  which  is 
larger  than  in  Figures  5.4a  and  5.4b  where  the  maximum  pressure  is  102  kPa. 

The  streamlines  are  superimposed  on  the  axial  velocity  contours  to  visualize  the  flow,  as  seen  in 
Figure  5.5.  It  is  clear  from  the  streamlines  that  none  of  the  exterior  flow  in  the  baseline  configuration 
(Figure  5.5a)  is  drawn  into  the  duct.  The  nozzle  configuration  results  show  that  the  velocity  is 
increased  downstream  of  exit  and  even  more  exterior  flow  is  being  drawn  into  the  exhaust  due  to 
entrainment.  Hence,  in  both  nozzle  configurations  (Figure  5.5b  and  5.5c)  as  well  as  the  baseline, 
the  streamlines  indicate  that  there  is  no  appreciable  flow  from  the  farfield  being  drawn  upstream 
into  the  duct. 
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Baseline  Configuration  5  Configuration  10  Configuration 

Figure  5.3:  The  duct  exit  of  the  three  configurations  in  this  study.  The  baseline  configuration 
does  not  have  a  nozzle  at  the  exit  (left).  The  nozzle  configurations  have  a  5°  (middle)  and 
10°  (right)  nozzle  at  the  duct  exit. 


Figure  5.4:  Contour  plots  of  the  static  pressure  (left)  and  vorticity  (right)  are  plotted  for 
the  three  computational  domains  analyzed.  The  baseline  configuration,  5°  nozzle,  and  10° 
nozzle  are  shown  on  the  top,  middle,  and  bottom,  respectively. 


55 


U,  m/s 


Figure  5.5:  Contour  plots  of  axial  velocity  with  super-imposed  streamlines  are  plotted  for 
the  three  computational  domains  analyzed.  The  baseline  configuration,  5°  nozzle,  and  10° 
nozzle  are  shown  on  the  top,  middle,  and  bottom,  respectively. 
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5.3  Computational  Grids  and  Boundary  Conditions 


Four  different  exit  boundary  types  are  considered  in  this  study.  Figure  5.6  shows  each  of  the  ge¬ 
ometries.  The  four  configurations  implemented  in  this  section  are:  imposed  back  pressure,  exhaust 
plenum,  and  two  experimental  configurations.  Experimental  configuration  one  represents  an  early 
experimental  geometry  of  the  rig  while  experimental  configuration  two  matches  the  current  exper¬ 
imental  configuration,  as  seen  in  Figure  5.1.  The  imposed  back  pressure  grid  contains  31,050  cells 
and  is  the  least  computationally  expensive  of  the  four  meshes.  This  configuration  models  only  the 
duct  and  does  not  include  the  exhaust  section.  The  boundary  conditions  are  mass  inflow  at  the  duct 
inlet,  no  slip  for  the  walls  and  flameholder,  and  outflow  at  the  exit  with  an  imposed  atmospheric 
pressure  of  101  kPa. 

The  exhaust  plenum  configuration  contains  58,342  cells.  The  height  of  the  plenum  is  twice  the 
length  of  the  duct  and  the  length  is  four  times  the  duct  length.  All  plenum  boundaries  use  an 
inflow/outflow  boundary  condition  that  allows  flow  to  either  enter  or  exit  the  domain  depending 
on  the  local  conditions.  A  pressure  gradient  is  imposed  across  the  length  of  the  plenum  in  order  to 
reduce  the  possibility  of  inflow  through  the  exit  plane.  The  pressure  on  the  left  is  105.5  kPa  while 
the  pressure  on  the  right  is  101  kPa.  The  cells  in  the  farfield  get  larger  away  from  the  duct  exit  in 
order  to  minimize  wave  reflection  off  the  artificial  boundary  surface.  This  configuration  is  the  most 
physically  desirable  domain  because  it  represents  the  flow  exhausting  into  a  large  plenum  section, 
avoiding  any  potential  for  downstream  influences. 

The  first  experimental  configuration  contains  50,828  cells.  The  farfield  has  a  17.78  cm  region  of 
inflow/outflow  with  a  pressure  of  105.5  kPa  immediately  downstream  of  the  duct  exit  followed  by 
a  106.68  cm  slip  wall.  The  wall  is  modeling  the  exhaust  duct  that  collects  the  flow  downstream  of 
the  experiment.  The  flow  is  exhaust  from  the  domain  through  an  inflow/outflow  boundary  with  a 
pressure  of  101  kPa.  The  pressure  gradient  is  added  to  reduce  numerical  instabilities  by  constantly 
exhausting  the  flow  downstream. 

The  second  experimental  configuration  contains  46,319  cells.  This  setup  has  the  same  boundary 
conditions  as  the  first  experimental  configuration  with  the  exception  of  moving  the  exhaust  duct 
in  the  farfield  to  be  60.96  cm  downstream  of  the  duct  exit.  The  wall  length  is  also  shortened  to  be 
39.37  cm. 
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IMPOSED  BACK  PRESSURE  EXPERIMENTAL  CONFIGURATION  1 


EXHAUST  PLENUM  EXPERIMENTAL  CONFIGURATION  2 


Figure  5.6:  Exhaust  boundary  conditions  and  computational  geometries. 


5.4  Effect  of  the  exit  boundary  condition  for  non-reacting  flow 

The  effect  that  boundary  conditions  and  domain  geometry  have  on  a  non-reacting  flow  field  was 
tested  on  the  rig.  Simulations  were  run  using  the  imposed  back  pressure,  the  exhaust  plenum,  and 
the  first  experimental  configuration  in  order  to  compare  how  different  outlet  boundary  conditions 
affect  the  results.  The  simulations  were  two-dimensional  non-reacting  flow  and  were  carried  out 
using  GEMS.  The  flow  was  driven  by  injecting  hot  air  at  2200  K  into  the  duct.  The  farfield  was 
initialized  with  atmospheric  conditions,  air  at  a  pressure  of  101  kPa  and  a  temperature  of  288  K. 
The  simulations  were  run  for  three  flow-through  times  so  that  the  results  could  be  analyzed  after 
the  initial  transients  dissipated. 

A  comparison  of  the  vorticity  contours  shows  that  all  three  cases  exhibit  Karman  vortex  shed¬ 
ding  off  of  the  flameholder,  as  shown  in  Figure  5.7.  The  pressure  contours  also  show  excellent 
correspondence  to  vorticity  for  the  three  configurations,  and  few  differences  among  the  configura- 
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tions.  It  appears  that,  at  least  qualitatively,  the  presence  of  the  exhaust  duct  or  the  imposition  of 
the  back  pressure  condition  did  not  affect  the  flow  features  in  a  significant  way  for  non-reacting 
flows. 

To  further  compare  the  three  domains  the  time  history  of  the  unsteady  pressure  was  recorded 
at  a  location  0.38  cm  downstream  of  the  flameholder  1.78  cm  off  the  center  line.  To  examine  the 
frequency  content  of  the  signal  a  power  spectral  density  (PSD)  analysis  was  performed.  A  total 
of  150  ms  of  temporal  data  were  analyzed,  with  the  data  spaced  1  fxs  apart  which  resulted  in  a 
maximum  frequency  of  2000  kHz  and  a  frequency  resolution  of  6.67  Hz.  It  can  be  seen  from  Figure 
5.8  that  the  spectral  content  for  all  three  configurations  have  strong,  well  defined  peaks  at  725  Hz, 
1450  Hz,  2175  Hz,  and  2900  Hz.  The  peaks  in  all  three  configurations  have  a  similar  amplitude. 
These  frequencies  correspond  to  the  Karman  vortex  shedding  and  its  higher  harmonics.  The  exhaust 
plenum  and  the  experimental  configuration  also  have  well  defined  peaks  at  600  Hz  and  1300  Hz. 
These  two  frequencies  correspond  to  the  vortex  shedding  at  the  exit  plane  of  the  duct  and  are  not 
found  in  the  imposed  back  pressure  case  because  the  artificial  boundary  condition  precludes  vortex 
shedding  from  taking  place.  With  the  exception  of  the  modes  due  to  the  vortex  shedding  at  the 
exit  plane  the  three  configurations  have  similar  behavior.  For  the  non-reacting  case,  these  results 
therefore  indicate  that:  (1)  the  experimental  configuration  reasonably  represents  a  large  exhaust 
plenum,  and  (2)  the  imposition  of  the  exit  boundary  condition  is  a  reasonable  computational  model 
of  the  experiment.  In  the  next  section,  we  examine  whether  these  conclusions  hold  for  the  reacting 
case  wherein  the  dynamics  is  not  caused  by  vortex  shedding,  but  by  thermo-acoustic  coupling. 
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Figure  5.7:  Contour  plots  of  the  static  pressure  (left)  and  vorticity  (right)  are  plotted  for  the 
three  computational  domains  analyzed.  The  imposed  back  pressure,  exhaust  plenum,  and 
experimental  configuration  one  are  shown  on  the  top,  middle,  and  bottom,  respectively. 
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Figure  5.8:  PSD  of  the  static  pressure  in  non-reacting  flow  at  a  point  located  0.38  cm 

downstream  of  the  flameholder  and  1.78  cm  off  the  center  line. 
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5.5  Effect  of  the  exit  boundary  condition  for  reacting  flow 


In  premixed  combustion  behind  a  flameholder  the  combustion  will  usually  suppress  the  von-Karman 
vortex  shedding.  Therefore  it  is  expected  that  the  spectral  content  of  the  pressure  signals  may 
be  very  different  for  reacting  flow  compared  with  its  non-reacting  counterpart.  The  prior  three 
domains  that  were  considered  for  the  non-reacting  study  are  repeated  here  using  a  reacting  flow 
setup.  In  addition,  the  second  experimental  configuration  which  has  the  exhaust  duct  wall  moved 
further  downstream  was  also  tested.  The  purpose  of  this  modified  configuration  was  to  further 
minimize  any  effects  that  the  exhaust  duct  may  have  on  the  reacting  flowfield  within  the  main 
duct.  Premixed  propane  and  air  were  introduced  into  the  domain  with  an  equivalence  ratio  of 
unity  and  a  temperature  of  288  K.  The  mass  flow  was  0.297  kg/s.  In  order  to  ignite  the  flow 
the  domain  was  initially  filled  with  warm  combustion  products  at  a  temperature  of  1792  K.  The 
simulations  were  run  with  the  GEMS  CFD  solver  using  a  single  step  global  reaction  [46].  The 
simulations  were  run  for  three  flow-through  times  so  that  the  results  could  be  analyzed  after  the 
initial  transients  dissipated. 

A  comparison  of  the  temperature  contours  shows  that  the  four  computational  domains  produce 
very  different  results,  which  can  be  seen  in  Figure  5.9.  The  imposed  back  pressure  domain  (Figure 
5.9a)  produces  a  more  uniform  temperature  profile  downstream  of  the  flameholder  perhaps  because 
the  artificial  boundary  condition  enforces  a  constant  pressure  at  the  duct  exit.  The  exhaust  plenum 
domain  (Figure  5.9b)  removes  the  constant  pressure  restriction  at  the  exit  which  is  the  ideal  sit¬ 
uation  since  it  allows  the  combustor  to  experience  farfield  pressure  effects.  This  result  shows  the 
presence  of  pressure  oscillations  and  a  pulse-like  appearance  in  the  temperature  contours.  The  orig¬ 
inal  experimental  configuration  (Figure  5.9c)  is  greatly  affected  by  the  pressure  reflections  off  of 
the  exhaust  duct  wall  in  the  farfield.  In  this  configuration  the  wall  is  located  17.78  cm  downstream 
of  the  exit  and  the  oscillatory  pressure  reflections  cause  the  hot  gas  to  travel  upstream  at  certain 
times  in  the  cycle.  In  the  modified  experimental  configuration  (Figure  5.9d)  the  exhaust  duct  wall 
is  moved  further  downstream  which  causes  the  flow  to  be  less  affected  by  the  pressure  reflections 
than  in  the  original  experimental  configuration,  however,  it  still  shows  some  differences  compared 
to  the  exhaust  plenum  case. 

It  is  clear  in  the  reacting  flow  cases  that  the  farfield  boundary  conditions  greatly  affect  the 
flow  features  in  the  duct.  To  get  a  qualitative  idea  of  this  the  pressure  time-histories  were  taken 
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at  the  same  point  as  in  section  5.4,  directly  behind  the  flameholder,  and  are  shown  in  Figure 
5.10.  The  results  again  show  that  there  is  a  significant  difference  among  the  four  domains.  A 
noticeable  difference  between  the  plots  is  the  amplitude  of  the  oscillations.  In  the  first  experimental 
configuration  (Figure  5.10c)  the  pressure  oscillates  about  20%  from  the  mean  values  compared  to 
5%  for  the  imposed  back  pressure  (Figure  5.10a)  and  10%  for  the  exhaust  plenum  and  the  second 
experimental  configuration  (Figures  5.10b  and  5.10d).  In  fact,  these  results  confirm  that  the  second 
experimental  configuration  indeed  represents  the  ideal  plenum  geometry  reasonably  well. 

The  PSD  analysis  of  the  pressure  time  history  is  shown  in  Figures  5.11  and  5.12.  A  total  of 
200  ms  of  temporal  data  was  analyzed.  The  data  were  spaced  0.5  /js  apart  which  resulted  in  a 
maximum  frequency  of  4000  kHz  and  a  frequency  resolution  of  5.0  Hz.  It  can  be  seen  that  all 
four  configurations  have  peaks  at  the  same  frequencies  but  with  different  amplitudes.  The  von- 
Karman  vortex  shedding  that  dominated  the  non-reacting  case  is  suppressed  in  the  reacting  case, 
as  expected,  and  the  acoustic  pressure  modes  are  dominant.  From  Figure  5.11  it  can  be  seen  that 
the  imposed  back  pressure  case  has  the  lowest  amplitude  which  directly  corresponds  to  the  smallest 
oscillation  from  the  mean  seen  in  the  time-history.  The  exhaust  plenum  and  the  second  experimental 
configuration  PSDs  agree  very  well  despite  having  somewhat  different  temperature  profiles.  The 
first  experimental  configuration  has  the  highest  magnitude  due  to  the  large  pressure  oscillations 
seen  in  the  time-history  plots.  Additionally,  Figure  5.12  shows  the  imposed  back  pressure  case  has 
a  high  PSD  amplitude  at  3200  Hz.  This  amplitude  appears  to  be  a  computational  artifact  that 
occurs  because  of  the  artificial  boundary  condition  imposed  on  the  domain  and  is  not  evident  in 
any  of  the  other  three  cases  that  involve  the  exhaust  plenum  section. 
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Figure  5.9:  Contour  plots  of  the  temperature  are  plotted  for  the  four  computational  domains 
analyzed.  The  imposed  back  pressure,  exhaust  plenum,  experimental  configuration  one,  and 
experimental  configuration  two  are  shown  from  top  to  bottom. 
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Figure  5.10:  The  pressure-time  history  is  plotted  for  the  four  computational  domains  an¬ 
alyzed.  The  imposed  back  pressure,  exhaust  plenum,  experimental  configuration  one,  and 
experimental  configuration  two  are  shown  from  top  to  bottom. 
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Figure  5.11:  PSD  of  the  static  pressure  in  reacting  flow  at  a  point  located  0.38  cm  downstream 
of  the  flameholder  and  1.78  cm  off  the  center  line.  (0-3000  Hz). 
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Figure  5.12:  PSD  of  the  static  pressure  in  reacting  flow  at  a  point  located  0.38  cm  downstream 
of  the  flameholder  and  1.78  cm  off  the  center  line.  (0-10000  Hz). 
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5.6  Conclusions 

The  type  of  exit  boundary  condition  was  examined  using  both  non-reacting  and  reacting  flows  in 
order  to  better  understand  the  influence  of  the  boundary  condition  on  the  spectral  content  of  the 
pressure  waves  inside  the  combustor.  The  study  was  also  used  to  help  design  the  experimental  rig 
by  evaluating  the  facility  effects.  In  particular  a  nozzle  was  added  in  an  effort  to  eliminate  reverse 
flow  into  the  combustor  from  the  ambient.  This  design  change  was  determined  to  be  unnecessary 
because  the  baseline  configuration  without  the  nozzle  showed  little  or  no  reverse  flow  into  the 
combustor  from  the  facility. 

The  non-reacting  flow  study  compared  the  vorticity  and  PSD  for  the  ideal  computational  (im¬ 
posed  back  pressure),  ideal  physical  (exhaust  plenum),  and  experimental  configurations.  All  three 
cases  produced  similar  results  with  the  only  difference  being  the  shedding  off  of  the  nozzle  exit  in 
the  latter  two  cases  which  does  not  seem  to  materially  affect  the  flow  in  the  reacting  duct. 

The  four  reacting  results  were  very  different  from  the  non-reacting  results.  The  temperature 
profile,  pressure  time-history,  and  PSD  analysis  are  heavily  dependent  on  the  outlet  boundary  con¬ 
dition.  The  presence  of  the  exhaust  plenum  (the  ideal  instance)  shows  that  the  temperature  contour 
and  pressure  oscillations  are  strongly  affected  by  the  external  flow.  Adding  the  exhaust  duct  wall  as 
in  the  experimental  configuration  creates  pressure  reflections  which  further  affects  the  solution  in 
the  combustor,  however,  it  is  less  dramatic  when  the  wall  is  moved  further  downstream.  Indeed,  the 
exhaust  plenum  and  the  second  experimental  configuration  had  similar  pressure  oscillation  mag¬ 
nitudes  and  similar  PSD  profiles  despite  having  different  temperature  profiles.  It  is  worth  noting 
that  all  four  configurations  showed  excitation  at  the  same  frequencies  but  had  different  amplitudes 
which  suggests  that  the  plenum  region  influences  the  thermo-acoustic  oscillations  in  the  duct.  It 
can  therefore  be  determined  that  it  is  necessary  include  the  proper  outlet  condition  when  modeling 
experiments  because  the  artificial  boundary  condition  (imposed  back  pressure)  does  not  capture 
the  correct  pressure  effects  and  also  shows  what  appears  to  be  non-physical  high  frequency  content. 
This  study  shows  that  boundary  conditions  and  facility  configuration  can  affect  the  computational 
solution  and  are  important  to  consider  when  performing  detailed  simulations. 
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CHAPTER  6 


Conclusion  and  Future  Work 


6.1  Summary 

This  work  analyzed  the  three  main  aspects  of  bluff-body  stabilized  flames:  stationary  combus¬ 
tion,  lean  blow-out,  and  thermo-acoustic  instabilities.  Stationary  combustion  and  lean  blow-out 
prediction  were  studied  using  an  improved  version  of  the  LEM  sub-grid  model  coupled  with  the 
GEMS  solver.  The  LEM  model  in  this  work  is  used  to  close  the  filtered  species  production  term 
in  the  large  scale  species  mass  fraction  equations.  In  previous  work  the  LEM  model  accounted  for 
the  species  mass  fraction  within  the  domain  solely  in  the  sub-grid,  where  as  in  this  version  the 
species  equations  are  solved  at  the  resolved  scale.  The  advantage  of  this  method  is  that  advection 
and  diffusion  between  LES  cells  is  correctly  represented.  The  results  showed  that  LEM-LES  over 
predicts  both  the  temperature  and  velocity  profiles  downstream  of  the  flameholder,  however,  the 
simulation  does  better  than  a  comparable  laminar  closure  model  in  capturing  the  temperature  field. 
The  LEM-LES  solver  did  not  predict  lean  blow-out  at  the  proper  equivalence  ratio.  Experiments 
have  shown  a  premixed  propane  flow  in  the  Volvo  rig  to  blow  out  at  an  equivalence  ratio  of  0.55 
and  in  this  work  an  equivalence  ratio  of  0.50  with  the  LEM-LES  solver  did  not  produce  blow-out. 
Therefore,  it  can  be  concluded  that  the  present  implementation  of  LEM-LES  solver  provides  some 
improvement  for  reactive  flow  predictions  of  velocity  and  temperature,  however,  the  model  is  not 
sufficient  to  predict  flame  extinction  correctly.  More  work  is  needed  to  improve  the  LEM  model 
in  order  to  better  predict  the  complex  dynamics  involved  in  turbulent  combustion  including  the 
exploration  of  more  sophisticated  chemical  kinetics  mechanisms  than  the  global  chemistry  used  in 
this  study. 

The  thermo-acoustic  instability  analysis  demonstrates  that  boundary  condition  effects  are  im¬ 
portant  to  considered  when  performing  both  non-reacting  and  reacting  simulations.  In  reacting 
flows  especially,  imposing  a  non-physical  boundary  condition  at  the  duct  exit  suppresses  the  acous- 
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tic  content  inside  the  domain  and  alters  the  temperature  contour  downstream  of  the  flameholder. 
The  addition  of  a  wall  in  the  exhaust  plenum  creates  pressure  reflections  that  also  effect  the  thermo¬ 
acoustics  inside  the  duct  and  is  necessary  to  consider  when  performing  simulations.  Therefore,  it 
can  be  concluded  that  it  is  important  to  include  proper  farfield  geometry  when  performing  calcu¬ 
lations  so  that  non-physical  boundary  conditions  do  not  influence  the  solution. 

6.2  Advantages  and  Disadvantages  of  LEM-LES 

The  LEM  model  is  attractive  because  it  has  validity  for  all  flame  regimes  and  it  solves  the  full  species 
equations  on  the  sub-grid  level  of  an  LES  grid.  The  model  is  able  to  capture  species  production 
and  diffusion  while  providing  an  reasonably  accurate  representation  of  turbulence  through  the 
triplet  mapping  process.  The  equations  are  solved  in  one  dimension  making  the  solution  much  less 
computationally  expensive  than  running  a  DNS  and  is  therefore  applicable  to  flows  with  higher 
Reynolds  numbers. 

The  LEM  model  is  computationally  expensive  because  it  solves  an  ODE  at  each  LEM  node 
within  the  sub-grid.  In  areas  where  there  are  large  gradients  there  can  be  as  many  as  thirty-six 
ODEs  being  solved  for  a  single  LEM  line.  It  is  estimated  to  be  about  five  times  slower  than 
a  comparable  laminar  closure  model  although  this  number  is  likely  to  be  problem  dependent. 
Additionally,  the  constant  pressure  assumption  in  the  LEM  model  makes  it  applicable  only  in  flows 
where  compressibility  effects  are  not  important.  This  limits  its  effectiveness  in  many  aerospace 
flows  and  acoustic  analysis. 

6.3  Future  Research 

Despite  its  advantages,  the  LEM  model  did  not  properly  capture  the  experimental  conditions  for 
lean  blow-out.  One  possible  reason  for  this  discrepancy  could  be  in  the  linear  initialization  process 
of  each  LEM  line.  After  the  LEM  solution  has  been  advanced  to  the  next  physical  time  step  the 
profile  is  discarded  and  initialized  as  a  new  linear  profile  during  the  next  time  step  (Reference 
Section  2.4.2).  This  does  not  allow  the  stirring  to  produce  a  proper  turbulent  profile  in  the  sub¬ 
grid.  In  future  research  the  LEM  profile  can  be  retained  after  each  time  step  so  that  the  turbulence 
from  previous  time  steps  is  retained  in  the  sub-grid  making  the  model  more  representative  of  a 


physical  flow.  In  doing  so  it  will  be  important  to  make  the  total  concentration  of  species  within 
each  LEM  sub-grid  line  match  the  resolved  species  concentrations  from  the  LES  solver  so  that  the 
model  is  physically  representative. 

The  LEM  implementation  in  this  work  treats  each  LEM  line  as  a  self  contained  entity  in  which 
there  is  no  advection  or  diffusion  between  LEM  lines  occur.  This  approximation  may  also  be 
contributing  to  the  erroneous  lean  blow-out  prediction  seen  in  the  Volvo  rig  geometry.  It  would 
therefore  be  important  to  further  the  model  so  that  these  physical  phenomena  are  properly  captured 
in  the  sub-grid.  Additionally,  future  simulations  should  be  run  with  a  finer  mesh  downstream  of  the 
flameholder  because  the  4mm  spacing  may  be  too  coarse  to  properly  capture  the  small  length  scales 
involved  in  this  problem.  It  would  also  be  interesting  to  run  the  simulation  in  three  dimensions  so 
that  the  flow  properties  are  better  represented  than  in  the  two  dimensional  simulations  presented 
in  this  work. 

Finally,  the  reaction  kinetics  used  in  the  Volvo  rig  simulations  was  a  single  step  global  propane 
reaction.  This  is  a  significant  simplification  because  it  does  not  account  for  the  many  possible 
associated  reactions  (perhaps  in  the  hundreds)  that  occur  in  a  physical  flame.  Future  work  should 
run  the  same  simulation  using  a  more  complex  mechanism  that  takes  into  account  intermediate 
species  such  as  radical  species  that  are  produced,  many  of  which  are  associated  with  ignition  and 
extinction  processes.  Incorporating  complex  reaction  kinetics  could  results  in  predictions  of  lower 
flame  temperature  and  therefore  more  accurately  predict  lean  blow-out. 
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